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ABSTRACT 

The semi*-'trayi*-by«»tray method for multistage, multi- 
component distillation column calculations provides the 
advantage of solving a variety of design and simulation 
problems with different end specifications to the desired 
degree of accuracy. The present work aims at overcoming the 
short comings of the available semi— tray— by—tray method by 
combining the tray-by-tray model of Naphtali-Sandholm and 
the shortcut model of Edmister with a simultaneous equation 
solving procedure. To promote convergence, the single variable 
Golden Section Opjtimizat ion technique has been coupled with the 
N&wton-Raphson method of solving simultaneous nonlinear 
algebraic equations. The linear flew profile approximation in 
the earlier semi— tray— by— tray method has been removed and the 
requisite flowrates are obtained by solving enthalpy balance 
equations along with equilibrium relations. Also, a more 
rigorous double stripping factor Edmister equation has been 
used as derived from the system of equations of the multitray 
model. The number of dependent variables has been reduced to 
a minimum still retaining the block tridiagonal structure in. 
the matrix which permits the use of efficient block Thomas 
algorithm. The generalized algorithm can be used both for 
design as well as rating of simple or complex columns. The 
progrsiii provides a choice of condenser types, tower configuration 
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and plate efficiencies. Ideal as well as nonideal thermo- 
dynamic property evaluation procedures have been incorporated 
in "the in case of nonideal mixtures^ the user has 

a c±iolce between SRK and Peng— Robinson methods or he can 
provide any other method of his choice. It also incorporates 
flexibility with respect to a wide range of end specifications. 
The CPU time and menory requirements are almost independent of 
the number of trays , 



CHAPTeR 1 


INTRODUCTION 

Multicomponent distillation plays a very important 
role in the separation process technology. Two different 
types of problems arise the simulation or rating problem 
and the design problem. In the rating problems, the separation 
that can be attained with a given column configuration under 
given conditions, is to be found whereas in the design 
problems, column configuration and operation to meet the 
desired separation specifications are to be determined. 

An efficient generalised method for solving multicomponent 
multistage distillation problems should essentially have the 
follcjwing features - (i) Simplicity of formulation, (ii) flexi- 
bility of application - i.e,, capability to handle problems 
with various complex tower configurations, feed conditions 
and locations, variety of end specifications and various types of 
eCndericeEii . and reboilers taking into account nonideality with 
respect to thermodynamics and stage efficiencies, (iii) Sound 
numerical technique ensuring stable and fast convergence having 
less memory and CFO time requirement, (iv) good structure 
and conditioning of the Jacobian matrix and (v) capability to 
handle both the design and simulation problems with desired 
degree of accuracy. 



Variety of approaches have been suggested for handling 
multicomponent, multistage distillation calculations such as 
short cut, tray-by—tray and semi-tray-by-tray methods. In 
the short cut method (Underwood 1932; Fenske 1932; Edmister 
1943 etc,), empirical correlations or a set of assumptions 
are used to solve the problem with comparatively less efforts* 
In the tray-by-tray methods either a sequential calculation 
is undertaken or the whole set of exact system equations is 
solved Simultaneously by an efficient numerical algorithm, 

Hiis is a rigorous method of solution and consumes larger CPU 
time and memory. Some of the important methods in this class 
have been discussed in a companion report (Das, 1985). The 
semi-tray-by-tray approach is a novel concept given by Ohm-.-n 
and Kasahara (1978), It maximizes the merits of both the 
above mentioned types of methods and still reasonably retains 
the accuracy of the tray-by-tray methods. 

The semi— tray-by—tray method, originally given by 
Olimura and Kasahara, had certain shortcomings and to overccme 

I 

tiiese, in the present study the semi— tray-by—tray procedure 
has been modified to include the tray-by—tray model of Naphtali 
and Sandho3jn and the multitray model of Edmister and a 
simultaneous equation solving approach has been used. The 
detailed formulation of the present method has been discussed 
in Chapter 3, Its applications to various end specifications 



have been discussed in a companion report (Das, 1985), In 
Chapter 4 the inccrporation of non*-ideal thermodynamics has 
been discussed. Some nearly-tideal and highly-nonideal systems 
have been studied and the results have been cited in Chapter 5 
A brief description of the nonideal subroutines in the 
computer package has also been provided. 



CHAPTER 2 


LITERATURE REVIEW 

A variety of methods have been reported in the 
literature to obtain the rigorous solutions for the two 
classes of problems namely •• design and simulation. In this 
review some of the prominent methods for design vjill be 
outlined. The literature on simulation methods have been 
rev5.ewed elsewhere (Das^ 1985), 

2,1 Advances in Design : 

Before attempting to review the algorithms available 
in literature to handle the design probloti, it should be 
made explicit that most of the so called "design methods” 
do not really conform to design and can only be economically 
and favourably used when the number of equilibrium stages in 
the tower and the feed locations are known. Exact solution 
of the "design problem" was first systematically approached 
by Lewis and Matheson (1932), They described a stage to 
stage calculation procedure to iteratively find the number 
of trays in the enriching and stripping sections which would 
allow a desired distillate product. This method assumes a 
product distribution and proceeds to calculate the component 
mole fractions at each stage simultaneously from the top and 
bottom of the tower. The calculation is terminated when the 



feed mole fractions are matched. An enpirical equation is 
used to generate the new corrected product distribution for 
the nesct iteration. A blo«k diagram for this method has been 
provided in Appendix A, This method, hcwever/has a nuiriJoer of 
shortcomings, which severely restrict its general application. 

It is essentially restricted to single feed distillation 
problems and difficulties arise in problems with extremes of 
composition. In superfractionator type columns, the top and 
bottom products are extremely pure. Under such conditions the 
system parameters may become sensitive to minute change in 
product compositions. Arbitrary procedures should be set up 
to handle non*-distributed components since the concentration 
are sensitive to initialization. Numerically also the 
sequential tray*»by««tray calculations are unstable due to building 
up of truncation errors. The method is susceptible to divergence 
for larger problems with nonideal mixtures. The rate of 
convergence by this approach is much slower than the more 
recent equation solving approaches. 

For a long time no new significant approach was 
developed in the field of distillation column design. In 1974, 
Ricker and Grens (1974) presented a calculation procedure which 
they named "Design Successive Apprcaclmation", This is essentially 
solving the design problem by repeated application of simulation 
procoduros. For the efficiency of this method, the number of 



column configurations examined must be small and accurate 
estimates of composition, flow rates and temperature must be 
computed at the loeg inning of each column calculation which 
makes the best use of values of previous column configurations. 
In an inner loop the Naphtali-Sandholm algorithm is used to 
generate the product distribution for an assumed number of 
trays and in an outer loop the tray numbers and Ng are 
changed to meet the split specifications. The change in the 
total stage requirement (i,e., Ng + Ng) is obtained from the 
difference between the calculated and specified reflux rates 
through a linearization of the Erbar-»Maddox; (1961) correlation. 
The corrected (Nj,/Ng) ratio is calculated frora a secant method 
convergence based upon the relationship between the number of 
stages in either section and the change in the key component 
ratio over that section as given by Pbnske (1932). Tlie 
remaining step involves the estimation of new temperature and 
component flow profiles for the altered column configuration, 
so as to preserve as much as possible, the convergence achieved 
in the inner loop in previous iterations. This involves 
scaling the individual ccmponent flows in proportion to the 
changed reflux and allowing for the change in the number of 
stages in a section, as well as scaling the temperatures in 
the middle portion of the column and holding temperatures 
unchan ged near the ends of each section, flow di art of this 

algorithm has been presented in Appendix B. Tliis method is 



reasonable for simple columns with regular profiles. However, 
attempts to solve the design problem through repeated simulation 
has run into several snags when applied to design of columns 
v/ith multiple-feeds or sidestreams or with unusual flow 
profiles, Ohmura and Kasahara (1978) proposed a new disti- 
llation calculation method utilizing the salient features of 
both the shortcut and tray-by-tray approach. This method can 
handle both design and simulation problems. Because of the 
novelty of the concept and since it also serves as the basis 
for the new proposed formulation, this method will be described 
separately in detail in the next section. 

2.2 Semi— Tray-bv-Tray Algorithm (Ohmura and Kasahara, 1978) : 

Naphtali and Sandholm (1971) method, while performs 
extremely well for operating problems, can not be used for 
design particularly design of complex columns with multiple 
feed and/or v;ith side streams, Semi-tray-by— tray method of 
Ohmura and Kasahara (1978) is claimed to overcome this 
difficulty. This method utilizes' both rigorous and short cut 
methods to obtain a solution close to trayWby-tray calculations. 
Hence the equilibrium stages are defined by two types of unit 
mod' Jls : 

A-Type Units ; where conventional tray-by-tray equations are 

used for each single equilibrium theoretical 


stage. 


f^nit s ; where shortcut calculation method of 

Ec3mister (1943) is applied and consists 
of an arbitrary number of trays in each group. 
The following assumptions are made for these . 
units, 

1, Loads (flot; rates) change linearly inside the 
unit , 

2, Single effective stripping factor by Edmister 
can be used in the system equation derived. 

3, Vapour and liquid leaving each tray are at 
bubble and dew point conditions. 

Figure 1 shove's how to represent a complex distillation column 
as an assembly of these two types of units. For theoretical 
tray model or Jwtype unit (shown in Figure 2) , the set of 
system equations are: 

1, Combined component material balance and 
equilibrium relations (Equation 6 in Table 1) 

2, Summation equation (Equation 4(b) in Table 1) 

3, Enthalpy balance equation (Equation 5 in Table 1) 
However^ it should be noted that in this case i refers 

to the unit number and not the plate number. It is also 
interesting to note that in a specific column with increasing 
number of A-type units, the calculation results approach those 
of the conventional tray-by»»tray method while with increasing 
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size of the B— type units the results approach those of the 
shortcut method. 


^♦2.1 Short Cut Model or the B«>tvpe Unit ; 

In the Birtype unit, M theoretical trays are dealt with 
as one model. Figure 3 shows a generalised Bwtype unit x^ith 
the variables. The variables with dashes indicate that 
they belong to the inside of the B«.type unit onlyw The system 
equations for the B^itype unit are 


^ ^ ^ Qysrall component material baiance for ^ the B'*»type Unit : 


^^1+1,3 + h-i,j - ''l,j <1+ 


(2.1) 


(b ) Overall Equilibrium Relation for the B-Type Unit ; 

Since in this method, an unit is being considered as a 
v/hole instead of a single tray,a relation betx/een the outgoing 
vapour stream (v^ j) and the outgoing liquid stream 
required. This relationship has been termed as the overall 
eciuilibrium relation for the unit 4 Nox'J', for any mth tray and 
jth component inside the ith B— type unit, the platewise equili- 
brium relation is given by 

V* .osS! jlj j (2,2) 

'^i/rn, j i/m,j 

K« ■■ V * 

where S! q « ■ (by definition) (2,3) 

^i,m 
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The overall equilibrium relation was deduced from 
tho component material balances and the platev;ise equilibrium 
relations by an approach similar to Edmister (1943)* Tliis 
deduction will be presented in Chapter 3, Edmister obtained 
an equation similar to equation 2*4 except an additional term 
which had a second stripping factor* Ohmura and iCasahara 
approximated Edmister* sexpression in terms of a single effective 
stripping factor without rigorously following the derivation from 
the relevant system equations , The final equation used by 
thorn is an approximate equation, using a single effective 
striiDping factor of Edmister, Se. ., given by 

1/ j 


i/j 


where D, ■* . 

j 


hv.J. 


(i+sb (1 -b D, ,) 
i 1/j 

se, , Cl - Se^'i*:^) 


(1^ W. 7) 

J 


(2.4) 


(2.5) 


and j = (Se, ,)'’U 

i/J i/j 


sM, 


Se 


i, j " C^i,l,j ^-25 J 


0*5 


With 


K* * . (1+S^) V 

S. = — A i 

1/1 


(2,6) 

0.5 (2.7) 

( 2 . 8 ) 


and 


s * 

i/M, j 


K ! V ! 

i*M,i i,M 

( 1+S^ ) iji^ 


(2.9) 


The difficulty arises when one proceeds to calculate 
^1,1 ^i^M* Ohmura and Kasahara made a major approximation 

that flow profiles inside the B^type unit are linear. Thus 
the vapour discharge rate from the lowest tray of the model 
v/as expressed as 



V 


( 1+S . ) -L . 

i+1 rnT^iT 



( 2 . 10 ) 


(c) Summation equation: (Equation 4(b) in Table 1) 

P*S!W Point and Bubble Point Equations : 

There are two temperatures in any B— type unit aiDpearing 
in the overall equilibrium relation through the respective 
equilibrium constants Kf - . and K! These top and bottom 

tray temperatures in the B-type unit have been defined by the 
respective dew point and bubble point equations as follovjs: 


1 

V. 


y 

j=l 


V 


iriL 


K) T 4 


0 


( 2 , 11 ) 


and 


c 

\ |i h,J - " ° 


(2a2) 


(e) Enthalpy equation ; (Equation 5 in Table 1) 


2,2,2 V ariables and Convergence Procedure ; 

The convergence procedure was similar to that of Tomich 
(1970) utilizing a two loop convergence scheme. In the first. 


loop, the component flow rates are generated by using the 
set of combined material balances and equilibrium equations. 

For B*^type units, a combination of equations 2,1 and 2*4 
have been used. In, the second loop, the enthalpy balance and 
the summation equations for the iwtype units (Equations 5 and 
4(b) in Table 1) were solved for and T^’s, However, for 

the B-»type units, the summation equation, enthalpy balance, 
defif point and bubble point equations were solved for top 
and bottom temperatures (T! . and T! ) and the numJoer of trays, 

X / X X / 

The convergence technique used was simple Newton— Raphson 
iteration v;ith direct inversion of the Jacobian matrix. In the 
first loop, however, a tridiagonal structure was present and 
heiiGG Thomas algorithm could be used. The correction vector 
x^as given by - . . 

■ £ j P (2.13) 

and the nev/ variables were computed as 

^NEW = ^ObD + (2.14) 

Ohmura and Kasahara have cited simple examples demonstrating 
the, application of this semi-tray— by-tray method, 

2,2,3 Shortcomings of Ohmura-Kasahara Method : 

Though this semi-tray-by— tray method is capable of 
handling simulation as well as design problems for even complex 
columns, it has several limitations. The method fails to handle 



complicated problems because of its list of simplifying 
assumptions and numerical instability. It^ thereby fails to 
satisfy some of the important criteria for a ;'eneralized efficicni 
algorithm discussed in Chapter 1. 

A critical analysis of the formulation shows that 
there is a large number of dependent variables/ specially/ all 
the comisonentwise vapour and liquid flowrates. The rigorous 
Jacobian calculation would involve finding the partial deri- 
vatives of these dependent variables with respect to the 
independent variables. Too many dependent variables lead to the 
difficulty of modulation inside the two loops specially for 
lariye problems. In a design problem neither the total number of ; 
trays nor the feed location are known. Hence the relative sizes 
of B— type units arc also unknown. In such a situation/ two | 

consequences are often possible, A convergence may not be | 

reached at all or if any shifting of dependent variaJales is | 

possible/ a distorted result may appear, A number of simplifying j 
approximations have been made which may not always be justified. [ 
A linearity of flow profile is objectionable for towers with I 

considerably large B— t;^’pe units as also in problems with i 

unusual profiles, Also/ in this algorithm/ provisions for ; 

nonideal mixtures and plate efficiencies x/ere not made. The ! 

authors have remained silent about the fle>cibility of the I 

algorithm and provision of end specifications to tackle various 


difficult problems. Numerically/ the algorithm seems to be 
poorly equipped. It lacks proper conditioning of matrisc and 
there is no sound matrix structure to enable the use of 
efficient numerical algorithms. The Newt on-Raphson technique/ 
wheii used by itself, can lead to divergence or oscillation 
for complex problems or unusual cases and therefore a guiding 
technique to convergence should be imposed to obtain a good 
convergence criterion. 

However, the algorithm has a very novel approach of _ 
incorporating ■ the number of trays as independent variables, 
Based on this concept, a new simultaneous semi*-tray-by'>»tray 
algorithm has been developed to overcome the difficulties 
discussed above and will be presented in the nect Chapter, 



TABLE 1 


Possible System Equations for Distillation Columns 
1, Material balance for jth component on ith tray: 


V. .. . + 1, \ V. .(1+s^) ^ 1. .(l-fsb + 0 


i+1 / j i«*l ^ j 1 / J 




irJ i/J 


2. Total material balance for ith tray; 

h+l + h-l - V^d+s^) - L^(l+s^) + = 0 

3, Equilibrium relation for jth component on ith tray: 


yi,j - ^i,j = ° h,j - h;j ''i,j = ° 


4, Summation equation for ith tray: 

c c - 

(b) y h;3' " h 

j=l j=l 


0 


5. Enthalpy balance equation for ith tray: 
c 


ti-i;j ^i-.i,j ^i+i^j ^i+i^j *“ ^i^'j ^v;j 

j=i 

*(l+sY) V, . K, . + . hf . + fY'h Hf'hJ =0 

i i,j 1,2 1,3 i,J i/j*'' 

6* Combined material equilibrium relation for jth component on 
ith tray: 

^i-l, j ^i+l,j ^i+l,j *" ^i,j ^i, j 


^i, j ■^-i/j 


o 
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7, Component material balance for ith tray including condenser 
rel5 oiler ; 


(a) -Aii/j. _ ^ 

%J b/j 


- 1 


(c) Note: for condenser (i,e,. i=l ) : = R + 1 

1/J 

f . 

(d) Note; for simple column: d, . ^ — 

l/J ■'-! 

1 +(^5^) 

^l,j 

Notati ons : 

■^i,J* -^^oEption factor for jth component on ith tray 

^i/j' P^o<^-tict dravjn from stripping section for jth 
component from ith tray 

d.. ' : product drawn from enriching section for jth 

component from ith tray, 

feed input for jth component on ith tray in 


'1/ J 


j 


P 


1* 


Xj v 

liquid condition (f. .) or in vapour state (f . .) 

1 / J ^ / J 

total feed on ith tray in liquid state or 

in vapour state (F . ) 


h^ ' j ; liquid phase enthalpy for jth component on ith tray 

H. j! vapour phase enthalpy for jth component on ith tray 

^i"j* pstt’tition coefficient for jth component on ith tray 

I, liquid flow rate for jth component from ith tray 

: r / j 
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; total liquid flow rate from ith tray 
: side stream fraction i.e., fraction of the stream 
taJcen as side jji^oduct, 

S. stripping factor for ith conponent on ith tray (= 1/A. .) 

3-/ J J 

v-"’.: vapour flow rate for jth component frcxn ith tray 
J 

; total vapour flow rate from ith tray 

Superscripts ; 

F . refers to the feed 

L : refers to the liquid phase 

V : refers to the vapour phase. 


CHAPTER 3 


hhaj semi^trj.y*.by~>tray method t^ith simultai^eous equation 

SOLVING (PRESENT METHOD) 

The ^erni— tray— by-. tray method o£ Ohmura and Kasiliara 
discussed in the previous chapter, includes a novel concept of 
combining the A and B type units xdnich lead to the possaiility 
of its use for design. The number of trays in a B-type unit 
can be _ unknown thus providing for its determination leading to 
design. However, there are some shortcomings discussed in the 
previous chapter, V7hich restricts its general application*. 

In the present method, Naphthali-Sandholm tyj,)e efficient 
solution procedure has been coupled with the smi-tray—by-tray 
approach. Also the assumptions e.g,, linear load variations 
inside the B«type unit etc, have been removed, We v;ill discuss 
these in detail. 

Figure 1 shows a multistage distillation column V7ith 
single tray units (A-type units of ohmura & Kasahara) and 
multitray units (B#*type units of Ohmura & Kasaii ara ) , The 
region with any disturbance such as feeds, side streams, condenser 
or rcboiler is provided with single tray units (models). The 
recjion of applicaloility of the single tray model depends 
entirely on the choice of the designer. Roughly one or two 



single tray units above and below the disturbance are recommended 
but it Is not a prerequisite. The remaining regions between 
the single trays are considered as multitray units, 

3,1 Single Tray Model ; 

It has already been mentioned that the single tray 
model and the relevant equations are similar to the ISfeiphtalaL-' 
Sandliolm method. It should, therefore^ be noted that i£ the 
number of trays and the feed location are known, as in the 
simulation problem, all the trays can be assigned single tray 
models and the computer package can be executed to obtain 
icigCMraus solution of the rating problem with Naphtaii-^andiiolm 
algorithm. Figure 2 rexoresents a generalised single tray 
model intth the symbols having their usual significance. 

The independent variables for each unit are: 

i) lj_^j j = 1,2,. ....,C 

ii) t. j = 1/2,, ,C 

J 

iii) T^ 

Here i refers to the unit number going from 1 to N and j 
refers to the component numbers going from 1 to C. There 
are a total of (2C+1) independent variables associated wfith 
each sincfle tray model* We need an equal number of equations 
to solve for the variables* These are: 



Ccmponent Material Balance Ecruations : 


1 “ -i ■“ 1. . - (1+sY) V. + V. , . 

i/j i—i/J 1 1/j 1 i/j i+l/J 


+ + f ^ 


where j=l/ 2 ,,.../C 


(3.1) 


^ Co mponentwise Equilibrium Relations with Murphree Plate 

Efficiencies, E. .: 

1 — 

E, , K, ■ . 1, . ■ V, , ■ ■ v . , T . 

— idl . -AU + (1-E. - .) (3.2) 


'i.j 


’^i ^i ^i+1 


where h 


1 ’ £ h., 

j=l 


. and V . = V . . 

3 i i/J 

j=l 


. and j=l,2/.../C 


and Kj_ j cite the equilibrium constants for the jth component 
on th e ith unit , 

(c) Entha lpy balance equations: 

G - 

E, = T fl. , . h. Y . + V. , . H, , - (1+sb 1. ,■ . 

i C-* L.i.»l,j i*-lyj 1+1, j i+l,j 1 i/J i/J 

j=:l 


(3.3) 

v/nere h denotes the enthalpy for liquid and H is the enthalpy 
for vapour and Q| is the heat supplied or taJcen out from 
each singletray unit. The sign of will be positive for 
heat supply to the unit and negative for heat supply to the 


unit and negative for heat withdrawal from the unit. There 
are {2C+1) equations associated with each single tray model. 

For condensers and reboilers these equations are 
adjusted according to the type of problem and given end 
specifications. The matching of the degrees of' freedom with 
the specifications and modification of the suitaJole equations 
have been discussed in detail by Das (1985), It should be 
noted that the equations for the single tray model are written 
in terms of certain functions such as M, Q. . and E. . These 
are called the "discrepancy functions" and have certain values 
or discrepancies when the numerical iteration continues and the 
unknowns are not solved. These discrepencies, however, reduce 
to zero as convergence is reached i,e., when the solution is 
found. 

The number of dependent variables is a minimum in this 
model. The dependent variables include only the equilibrium 
constants (K. .) and the enthalpy terms (h. “". and H.- , ), which 

XfJ l/J 1,J 

arc dependent on temperature and pressure only for the ideal 
case. For the non-*ideal systems, comparatively more camplicated 
functionalities arise and those will be discussed in Chapter 4, 

VJhen all the N-trays in the 'tower are taken as single 
tray units, the formulation is essentially that of Naphtali.-' 
Sandholm type and all the N(2C+1) equations can' be solved 
tocjother to obtain the colunn performance. It should be noted 
here that all the (2C+1) equations for ith unit have only three 
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sets of variables i.e., variables associated witli (i~l)th 
unit/ ith unit itself and (i+l)th unit. This gives rise to the 
syaimetric tridiagonal structure in the matrix. We v/ill proceed 
to discuss _ the matrix structures after describing the multi- 
tray model. 


3,2 Multitray Model : 

.'.■■'igure 3 depicts the schematic representation of a 
general multitray unit along with the incoming, and outgoing 
energy and mass streams. In any ith multitray unit, as sha-ni 
in i'igure 3/ Mj^ trays are incorporated and this M^ can be 
treated/ as an independent variable for design computations. 

The independent variables, for , each , multitray unit are: 


(!) 1 
(ii) V 
(iii) M, 


i.j 

i.j 


J — 1/2/,,,, ../C 
j = 1/2/ ./C 


Mere i refers to the unit number going from 1 , to N and j 


refers to the component nurtber going, from 1 to C, There are 
a total of (2C+1) independent variables associated with each 
multitray model. We need an equal number of equations to solve 
the variables. 

The equations used in the multitray model are as 


follows : 



( ^ ) ^^JPonGD Mat.G]rjLal BalancG Equations * 

"l.j = - (l+sj;) - (l+s|) v.;j (3.4) 

Where J = 

(b) En thalpy- Balance Equations ; 

E. = y ll. T -■• . h, - •. -d+s?) 1. . h.". - (1+sY) V.". H. 

jImp ^ i*«'X X i/j X. x^j x^Y) 

j=i 


+ h+i,j + 4,i * 4 


- . hY- .T + Q? (3.5) 
/j t/jJ 1 


The material and energy balance equations are identical in 
case of multitray unit and single tray unit. However^ the 


cliffejccn af3 between the two models appears in the eejuilibrium 


equations v;hich follow;' 

(c) Overall Equilibrium Equation for the Multitray Unit : 

The overall equilibrium equation is obtained by 
condensing the platew.lse ccmponait material balances and 
equilibrium relationships and eliminating the varialoles inside 
the multitray unit as shoX'Jn in the following deduction. It is 
knotni that equilibrium relation and component mate3rial balance 
for a general mth tray in the multitray model are given by 


+ 1 ! 


= . V 


+ 1 ! 


i/ft/j i/m,J i,m+l,j i/m«.l^j 


V 


S‘ 


. 1 ', 


i,m^j i,m,j i,m,j 


(3.6) 

(3.7) 





VariaDles with primes indicate that they refer to trays 
inside the multitray units. 

Starting from one end of the multitray models have 


i/l.j 


1,1, j 1,1, j 


S' - (vi _ + 1 . - vl'v"-) 

1/1/1 1/2, j i-l/J 1/1,1 

(Using equation 3,6) 


1 ■! sro 4 1| o • + Sj . (1. - . v! , .) 

1^1 ^ i ^ 1 ^ j 


(Using equation 3,7) 




e I . . . . 11.., . , . 

1,M,J 


(S' - . + S! - . S' ^ ^ 

i,l,J 1,1,1 1-/2, j 


+ 


.+ 


'l Sj o 4 

1/1/1 1/2,1 




Ijet 


-r. 


-s q I cl ^ Cl - . 

i,j ±,l,j i,2,j i/M,j 


(3,8) 

(3.9) 


and 


D. 


+ S' 


Sj n , i 


i,j “ i/l,j i,l,j i,2,j 


.+ 


®i,l,j ^i,2, j •••" •'•^i/KUl^ J 


(3.10) 


Also, it is obvious (Figure 3) that 


i,l, j 


(1 + sY ) V 


i' "i/ j 


(3.11) 


^i,M/j ~ ®i^ ^i,j 


and 


(3.12) 



Ohmura and Kasahara deduced equation 2.4 by substituting 
emations 3,9., 3,10, 3.11 and 3,12 in equation 3.8. In the 
present formulation a more rigorous double stripping factor 


approach of Edmister has been used and a modi.fied system 

equation has been formulated from which the values of the 

stripping factors cani be directly dedoced* 

From equations 3.8, 3,9 and 3,10 we have 

V ! ■ , . = .11 . + ( d . . ) (1 . - . - v ! - ■ .) 

1^1/1 i/j i/h/j i/j i-l/j i/l/j 




1 ! 




(D. + IT--,.) (1. 


: 


) 


or 


■' -■ (1+w, 






^^rhero (SUM) = * ^i-l,j ^i,l,j^ 

and 


or 


or 


W, .= (D. , + TT .) 

irl 1/J 


(3.13) 

(3.14) 

(3.15) 


W 


V 


iti 


1. 


i/l/j (i+w .) 

1 / J 


w. . 

1 , 


■iT'i-'. 

+ (SUM) 


(1+W. -‘ .) 

J 


(1+W, ,) 

J 


+ (SUM) ~ (SUM) 








V ! . 

i / 1 > . 1 


= (SUM) + 
-(SUM) 


1+W 7f , ,, 

1 - . - (SUM)l duAT 

I*"-!- / J 


w 


i/ j 


w. . 
1/ 1 
(T$^ 

a. 


^i*-l , j 


1 - 


(SUM) 


1 . , .Si. 
l-l/j 1^1 


w. 


(3.16) 


(1+I'h .) 
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Using eqn. 3.14 in eqn. 3.16 we have 


1. , . - 1> 
— 




(SUM) 


1. 1 • . S' . 
1-1 /J l/J 


w, 


i/JL 


( 14.W . . ) 


where 


S' . 
1/J 


VJ . . 


1+W.',«J7. . 

i/J 1/3 


(3.17) 

(3.18) 


Let Se. ‘ . represent the effective stripping factor for the 

X/ J 


jth component in the ith unit. 

M. 


Therefore ; 


and 


and 


and 


TT = se .'i , 

i/J i/j 

°l, j - T s^, , '■ 


- ®®i,j + ^®i,J + 


Mi 

■^®i/j 


M 


JL+l 


(Se. . - Se. j )/(Se . - .-1 ) 

X/J X/J 

M - 

+ + 

M.+l ■ 

(Se. . •*> 1 )/ ( Se . . •* 1 ) 

i/j i/J 


(3.19) 

(3.20) 


( 3 . 21 ) 


(3.22) 


Nq./ substituting equations 3.19 to 3.22 in eqns. 3,17 and 3*18 
v/e have ; 


^i-.l^j 


(SUM) 


^i«l,j ^i/j 


Se^^it^ - Se. '. 
l.> 3 : 

S„^ - 1 

i/j 


( 3 .23 ) 



Edmister (1943) in his deduction used the assurnption 
that most of the mass transfer takes place in the terminal 
IDlates of the column only and hence the expressions for the 
ef f ectivestripping factors can be deduced using a column 
with tv7o trays (i,e,^ by taking the first and Mth trays only) 


So. 


e. . 


- S. 


.1 


ll 


'i/ j 


Cl .. 4- ^ * q I u. 

i^l,j ^ i,l/J i,M/j 

(S‘ . + S'. - • . Si .+1) 

1,1, j i/l/J 


(3.24) 


By factoring out the term (Se. . -1) and simplifying we have 

J 

the quadratic equation 


S§f . + Se. . 

i/j i,j 




wliich has the positive solution 

~ t^±,l,j 0.25J 

Similarly, from equations 3.18 we have 


0.5 


0.5 (3.25) 


W. ■ . 

JuA. 


S! , ■ . + S! . SI . . 

-t/l / j 






* S ' m, .) 
1,M, j 


(3 .26 ) 


S j is the second stripping factor term of Edmister. 

sujDStitutirg equations: 3.11 and 3,12 in equation 3.23, 


we have: 





d+si') 1. . 





'l*-! j j 


1 - 


CSIM) 


1 . , . S ! - . 

1-1, j i,j 


where (SUlvi) = (l+s^')l. . « 1. ^ , + (1+sY) v. ' . 

1 i/J i-ld 1 i/J 


Mj+1 

Se . ' . - Se . . 

1 , j a. , j 

ir+i~= 

Se. . -1 
i/J 

(3.27) 

(3.28) 


This is the rigorous form of the overall equilibrium equation 
to be used for the multitray units in the convergence 
calculations. 

It should be notSd here that the expressions for the 

effective stripping factors consist of two terms SJ-. and 

1/1/1 

which are given by: 



Ha,! = 

(3,29) 

and 


(3.30) 


and here we have tv7o dej^endent variables L| ^ and Vj^ which 
have not been defined so far. These variaJoles can be associated 


in two ways, 

(1) Approximate m et hod : In this approach, the flow profile 
inside the multitray model is assumed to be linear, Tnis 
necessarily implies that the flowrates on successive trays 
increase or decrease by a constant value throughout the 
multitray model. That is, we have 



31 


L! ■, 

1,1 


V! 

i,M 


Li . ^ 

1-1 


4 (1 + 


M. 

1 * 

-v^(l+S^) - 

-M. 

1 


+ L 


'i-1 




(3.31) 

(3.32) 


For very simple problems, this method can be used to 
simplify computational efforts and thereby save time. However, 


only an approximate solution is obtained in this case and the 
convergence cannot be guaranteed. For complex, nonideal cases 
or problems with irregular flow profiles the rigorous approach, 
to be discussed, is suggested, 

Vigorous Method : In this method, and ^ are 

calculated rigorously from the system equations using mathematical 
techniques. 

The general dew point and bubble point equations for the 
Gomponen'- s leaving the ith unit are given by , 

V. . ‘ 

-iU ^1=0 = f(T, ) (3,33) 

JfN. , , 1 

i/J 



and 



V' 

HI I ■ 

j=l 




1 = 0 = f(T^) 


(3,34) 


Now, for the multitray model, considering the Murphree tray 
efficiencies to be unity, the vapour leaving from the top is 
at the dew point temperature, Tj^ ^ and the liquid leaving from 
the bottom is at its bubble point temperature, If the 

coraponent vapour and liquid rates leaving the multitray unit 
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axe known, be calculated from equations (3,33) 

and ?.34 r'^.c-cctively. Consequently, the mole fractions of 

the liquid leaving the first tray in the multitray model 

(i,e, 3cJ .) can be obtained from; 

/ J 



V. 

X 



(3.35) 


Tlio enthalpy balance for the first tray inside the multitray 
moc'«";l is given by 


G 

_ 




j=i 


vl - , H! - , + 1! T . h*. T j ■“ Ij T ■ j T ' 

1/1/1 1/1/ j i“l/j i*"!/! 


•v! „ . H! _ , 1 = 

1/ 1 1/^/1 '( 


0 


(3,36) 


Also 


II S= -Jf • L' 

h,i,j 1,1 


(3.37) 


and 


v!- .=v! T .+1! T .-I. . . 

1/2,1 1/1/1 1/1/1 1-1/1 


(3,38) 


Therefore, substituting equations 3,37 and 3,38 in equation 
3,36 Tfie have 


"S” f(l±sY) V. , H' - . + x* T . L‘ , hi', ■. - 1. 1 . h^^^ . 

b. - i i,j 1,1/1 1/1/1 1/1 1/1/1 1--1-/ J -T---L/J 

j=l 

- f j + ^1,1, j h,i - h-i,j t d ,2,j J “ 

' ' -/■ . 

(3.39) 


o 
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H’ , the enthalpy term for vapour leaving the second 
tray multitray unit, can be computed if the plate 

temperature ^ is known, which, in turn, can be calculated 


from an equation (Eqn. 3, 40) similar to Eqn, 3,33, since v! 


are known as functions of L! . , 

1, i 


1 

V 


• Z « 1 = 0 = f (T* - ) 

.J iry i fc wi m Jf ^ JL / « 


i / 2 , j 


(3,40) 


i.2 ‘'i,2,j 


Hence the nonlinear equation 3.39 can be written in the form 
f (Lj^^l) = 0 from which L| ^ can be numerically solved using 
Regula-falsi technique. Similarly, by using the component 
material balance, equilibrium arid enthalpy balance for the 
bottom plate of the multitray mode]^ M computed. 

The equations are as follows. 

Enthalpy balance for the bottom tray inside the multitray 
model is given by ' ' . 

c ' ' , . .... 




fji 4. 1 » ?i • ' 1 • • h *■ 

M/J i/^/j ■ i,M,j •^i,M-l,J i,M-l,j 


j=l 




Where 


V 


y\ 


v 


i/M,j ■'i,M^j i,M 


■ (3.41) 


(3.42) 


and 


"1 I — ( \T *4“ 1 ® '■* 




(3.43) 




with 


TL K! ,, 1,. . - 


: 


j=i 


j i/ j 


1=0= f(T';„) 


(3.44) 


and 


G 

■'ST 

K! 


. 11 




(3.45) 


3,3 Variables Selection and Matrix Structure; 

It has already been shown that each single tray unit has 
(2g+1) independent ttariables and equal number of equations. 
Retaining the consistency of size and symmetry in structure of 
the Jacobian matrix for all the units in the column is important 
for CPU time and memory reduction* This essentially demands 
that the number of independent variables be restricted to 
(2C+1.) only. The output variables, i.e*, vapour and liquid 
flow rates issuing out of the multitray units already constitute 
2C independent variables. 

Further, the multitray units have two terminal temperatures 
to be considered per unit in addition to the number of trays, 

If all these are treated as independent variables, the total 
number for multitray units will be 2C+3, This will destroy the 
constancy of size in the Jacobian matrix, the number of variables 
for a single tray unit being 2C+1, To overcome this problem, 
only the number of trays, has been taken to be an independent 
variable. The terminal temperatures along with the total 
vapour and liquid rates discussed earlier have been considered as 
•Tapendent variable®- and are computed, whenever requited, from. 
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the set of independent variables. To sum up, following have 
been selected as independent and dependent variables for 
multitray unit; 

Independent Variables 

(a) coraponent vapour flow rates 
from the top of the unit 

(b) Component liquid flew rates 


from the bottom of the unit 

(g) Nunaber of plates inside the 
unit 


Variable Name Number 








The number of independent variables is 2C+1 per unit which is 
equal to the number of independent descrepencies specified for 
the unit, 

B , Dependent Variables : 

(a). Total vapour flow rate from 
the top plate of the unit 



(b) Total liquid flow rate fran 
the bottom plate of the unit 


(c) Top plate temperature of the unit 


c 



j=l i 

T!"-, evaluatSkl from I 

JL ^ X t: 

the dew point equation 
of the outgoing 


vapour 
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(<a) Bottom plate temperature 
of the unit 


(e) Top plate total liquid 
flow rate 

(f) Bottom plate total vapour 
flow rate 


T! evaluated fron the 

bubble point equation 
of the outgoir^ liquid 


L! • evaluated as mentioned 
1^1 

in Section 3,2* 


V! evaluated as -mentioned 


in Section 3,2 

There are in total, N(2C-fl) independent variables and an 
equal number of equations where N is the total number" of units 
in the column. Now we proceed to introduce the vector notations 
for the column matrices. Consider as the vector of variables 
for unit i -and as vector of functions for the same unit. 


X = 




I 

HI 

4m 


mt 

^2 


^2 

•> 

and p = 

t 

1 

t 



1 


1 ^ 

t 

\ 



M 


N’ 





(3.46) 



o / 



J, the Jacobean matrix which is really a matrix of the partial 


derivatives of all the functions with respect to all the 
variables at the present value of the x ^ is given by 



( 3 . 48 ) 
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v;lth 




d ^j,2C+l 


^ ^j,2C+l 


(3.49) 


Fortunately^ the distillation column equations contain 
’ "iriables pertaining to three consecutive plates/ namely, (i“l),i 
and (i+1) only and therefore, in the above matrix most partial 
derivative terms are zero except three and hence the Jacobian 

! 

i 

matrix ai-sumes the block tridiagonal form as follows: 



(3.51) 
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Now we proceed to determine the structures of tJiQ 
Rubmatrices B\ and in the overall Jacobian matrix. The 
structure for the single tray model will be similar to the 
Naphtali^iSandholm algorithm and is shown in Figure 4, For 
the multitray models, however, there will be a change in the 
structural pattern as shown in Figure 5, The cross marks (X) 
in the pattern indicate the presence of the derivative term 
at that location. However, this pairtial derivative can be 
calculated either numerically, analytically or in a mixed mode 
so as to minimize computer memory and CPU time requirements, 

3,4 Solution Procedure: 

Neiwton*>J^aphson solution procedure for a set of nonli;: 
equations requires - 

« 5?^ + ^5?^ - (3.52) 

where I and I+l are iteration numbers and is the 

correction vector, given by 

= - 1^ ^ (3,53) 

In order to evaluate the correction -vectcr, one must 
compute inverse of the Jacobian matrix 5*^, Of the various 
procedures ®vailaible,Gaussian Elimination and Thomas Algorithm 
are the most important. Since the structure is block tridiaoonal 
the Block— Thomas Algorithm is the most efficiert « The user 
has been provided with the option to choose either of the two 
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procedures in the present prcjgram (Appendix. F)* A detailed 
discussion on these algorithms and their CPU time and memory 
saving capabilities has been provided in a companion report 
(Das^ 1985). 

However^ it is not guaranteed that the correction vector^ 
when used as such^ will ensure convergence. There can be a 
number of reasons Which can lead to a diverging or oscillating 
behaviour. For example^ the convergence is sensitive to 
initial guesses and the local behaviour of the Jacobian matrix. 
Therefore, it is suggested that a fraction of the correction 
vector is taken and equation 3,52 should be rewritten as 




( 3 . 54 ) 


Where y lies between —1 and 1 (Kaufmann, 1983). To implement 


equation 3, 54, it is absolutely necessary thst p is selected 
in such a way that convergence will be ensured. To achieve this 
criterion, a single dimensional Golden Section (Beveridge and 
Schechter, 1970) search routine has been incorporated* 

The unidin ens ion al numerical optimization f'Car unimodAl 
functions using Golden Section technique involve the following 
steps. 

Step 1 : Define the range of search for f i,e., ^ ^ 4 

Step 2s Locate two points ( and ^ 2 ) ^ distance of 

(l-l) and 1 from one end of the total range, 
n/s-I 


where 1 » 


* 0,6180 (approx,) 
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Step 3 ; Calculate tl>e variables X 

Step 4: Repeat Step 3 for 


I+l 


from equation; 3*48 


f 


Step 5 : Calculate the overall sum of the squares of all 


P * ^ 1 ^^[<3 ^ = 1 ^, 


discrepencies,, (SUM)j^ and (SUM)^ for 
respectively. , . 

Step 6: Calculate (SUM)^ and (SUM)^ corresponding to the 

jwT4-l 

va^lue of the variables a from equation: 3,48 with 




and 


f-h respectively. 

Step 7: A new interval of sea^cch .is selected, thus; 

(a) If (SUM)j^ equals (SUM )2 the ne*w interval for P 


becomes 


fl 


(b) If (SUM)j^ is greater than (S[JM) 2 , the new interval 
for ^ becomes 

(C) If (SXJM)^ is less than (SIM) 2 , the new interval 
for ^ be cone s 

Step 8: Define appropriately and ^ the boundaries of 




, according to the 


of the new range of search for 
particular case and repeat Step-2 to Step— 7 till the 
absolute value of ( 

where is the supplied tolerence limit 

It should be noted that for most of the ideal cases the 


3^ -^^ 2 ^ is less than or equal to ^ 






value of ID tends to unity. However, for complex situations 
assumes, different optimal values to guide the iterations towards 
convergence. 
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Figure 5: Structural pattern In the submatrices A.^B, 
and 6f the overall Jacobian matrix whel 
i represents a multitray unit nurrtaer only 








CHAPTER 4 


THERMODYNAMIC PROPERTY ESTIMATIONS FOR MULTICOMPONENT 
DISTILLATION CALCULATIONS 

4,1 Introduction i 

A prerequisite for any simulation (or design) 
calculation of a chemical plant or process equipment is the 
knowledge of physical and thermodynamic data of various streams 
(pure components and their mixtures) in the process. Almost 
all good simulator packages contain an internal library of 
thermo-physical data and also have inbuilt programmes for the 
estimation of thermodynamic and physical properties. 

To solve the distillation problems involving multi- 
component mixtures/ the vapour— liquid equilibrium data and 
enthalpy data are required. The methods used to obtain these 
data may be classified as follows; 

1, Curve fits of pure component data obtained 
experimentally or from generalized coixelations 

2, Generalized correlations as functions of temperature/ 
pressure and composition for both pure components 
and mixtures. 

3, Techniques using a single equation of state for 
both liquid and vapour phases. These are known 
as "Direct Methods", 


4, Procedures using an equation of state or 

ccrrc.lrtion for the vapour jAiase and a group 
contribution method for the liquid phase. These 
are known as "Indirect Methods", 

The first category methods are used when a large amount 
of experimental data is available. This is usually recommended 
for very odd cases having association and other interactions 
or for uncommon systems with unestablished thermodynamic 
behaviour. The methods of the second category utilise some 
correlations, which have been repeatedly tested and validated 
for a large number of systems. These are extremely useful for 
ideal mixtures, however, they often pose problems when dealing 
with nonideal systems. Expressions, vjith respect to temperature 
pressure and vapour and liquid phase mole fractions, having 
generalized applicability to all systems at all ranges of 
conditions, are very difficult to obtain. Examples of these 
methods are Antoine Equation (1888)? Yen and Alexander (1965) 
correlations etc. The third category of methods emerged 
recently when it was found that the same equation of state can 
be applied to both liquid and vapor phases in most of the common 
systems with fair degree of accuracy. However, the most 
commonly used equations of states sometimes fail to predict 
the liquid phase nonidealities of highly associated or inter- 
acting liquids. The industrial uses of the equations of state 
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Ui’ 

have been elaborated in paper by Adler et al,(1977). 

For complicated industrial systems, the use o£ multiple equations 
of state for liquid and vapor phases has been suggested by H 9 llan< 
(1981), These equations have been elaborated in Section 4.3, 

Fcsr systems having highly interacting or associated behaviour, 
a single equation of state often fails to predict the behaviour 
in both liquid and vapour phases. In such cases, the fourth 
category o-f methods i,e., the indirect methods is recommended. 

It utilizes molecular interaction or group interaction parameters 
to predict the liquid phase activity coefficients. This method 
is particularly useful to associated systems like alcohols 
and will be presented more elaborately in Section 4,3, However, 
it should be mentioned that the degree of accuracy required 
and the exact method to be used is upto the choice of the 
designer and a trade off between accuracy and computer timing 
should be mad^ iti accordance to the requirements, 

4,2, Basic Property Est i mation for Pure Components : 

The basic properties which are required as an information j 
feed to the rigorous nonideal calculations, include the physical 

■ I: 

constants like molecular weights, boiling points and critical | 

' ! 

properties for pure components. For very common compounds 
it is suggested that those values can either be directly given 
as an input or may be automatically read from a property data 
bank, HoxTOver, for uncommon systems like halos ubst it uted ! 


i: 



silanes very little information is available in literature* 

The reported experimental data for such unconmon systems has 
an uncertainty as high as 40% to 60%. For such cases estimation 
of properties from the molecular weight and boiling point 
data is essential* 


4,2,1 Critical Propecty Estimation : 

Critical temperature^ pressure and volume (i*e* 
and V^) represent three widely used pure componcait constants* 

The critical properties for most of the common compounds 
have been given by Reid et al. (1977), Comprehensive revievjs of 
critical properties are available by Kudchadkar et al, (1968) 
for organic compounds and' Mathews (1972 ) for inorganic compounds, 
I'Jhen estimation is required", the method of Lydersen (1955) is 
normally employed. In this method, structural contributions 
for temperature, pressure and volume (i.e, and Ay 

respectively) are employed to estimate T^, and V^. The 
relations are 




-1 


M (0,34 +y Ap)' 


(4,1) 

( 4 ’. 2) 


and = 


Uq (4, -3) 

where T. = normal boiling point of the compound in degrees 


and 


M 


Kelvin 

molecular" weight of the compound* 
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^ A A and for various compounds have 

Values of At# aod y 

al (1977) as well as by Lydersen (195 
b^en'i given et a . 

T , P and V, are obtained in degrees Kelvin, aton,ssberes 
and cli/centi^eters per gra„ .ole respectively, errors 

associated in estimation are usually less than TA for and 
lass than 4% for P,. Av <»^ta for all compounds are not 
available for the lydersen's m^hod. So for seme organic 
ooapounds. .. mo^iod of Riedel (19S4) is .uite accurate. In 

this method 


3,72 + 0.26 ( 


where 


= 0 19076 1 1.0 + 




1,0 - T. 


(4.4) 


(4.5) 




With P, in atmospheres. R isthe universal gas constant, 
82,06 cm^ atm/gmmole K 

OC is kncjwn as the "Riedel factor". 

- . 4, 4-vio critical compressibility 

At the critical pointy the cnri^aj- 


f act or, Z^, is given by 


c c 


(4.6) 


4 . • rac -for the halogen—substituted 

The estimated critical properties for th 

t- ,,v, in Table 2 as an example, 
silane system have been shown in Table 

4.2.2 Estimati on of Acentric Pagto r. 

one of the most common pure component, constants iS the 

/TD-n-^c-r 1955 ), ’Which is defined as 
Pitzer acentric: factor, ( ■ * 
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« log P^p ^ (at = 0,7) - l"o (4,7) 

where T^e '^c/'^b and ■ (4,8) 

^VP r” ^VP^^c ^^VP vapour pressure) (4.9) 

As originally proposed,. cs was to represent only the acentaricity 
nonsphericity of a molecule. For monatomic gases, is, 
therefore, essentially zero. For higher molecular weight 
yydr ocarb ons , increases and often rises with polarity. At 
present it is widely used as a parameter, which in seme manner, 
measures the complexity of a molecule with respect to both the 
geometry and polarity. It should be noted that for strongly 
polar or hydrogen bonded-fluids, the above correlation is 
liable to involve errors. 

To estimate the acentric factor by equation 4,7, the term 
Pyp is required, which, in turn, depends upon the vapour 
pressure equation used, 

Edmister (1958) suggested 


B 


log Pyp = (A + ^) where A and B are constants 


and 

xvith 




3 

7 


r g T 

11- ®J 


log P^ -1 


(4.10) 

(4.11) 


» = V^c 

Lee and Kesler (1975) suggested an improved correlation f or 


The relation is 



-In P « 5.92714 + 6.09648 + C28862_ln@ - 0,169347$^ 

SS 

15.2518 - 15,6875 ®“^- 13,4721 In© + 0.43577 O 

(4,12) 

The estimated acentric factors for the halosubstituted silines 
systems have been reported in Table 2, 

4,3 Calculation of Partition Coefficients, K,- . : 

- .... - ^^ 3 

In distillation calculations, the nonideality of thermo- 
dynamics enter mainly through the partition coefficients or 
equilibrium constants for the ith stage and jth componcnt,K^^'j, 
If the system to be studied is ideal or near ideal, the 
can be approximated to be a function of temperature T^^ and 
pressure P only. However, if the system is nonideal, the 
apart from the above functionality, becomes a function of the 

liquid and vapour molefractions (x, • . and y.* ■ .) also. For 

1# J J 

convenience sake, in this chapter, we will drop the subscript 1 
and denote the equilibrium constant as K^, Now we proceed 
to the estimation procedures-. 

4.3,1 Estimation of K.’s for Ideal Systans : 

■■ - I.'.* - .■WI . II W I ’ ■ ■— i. III 

As already mentioned, for ideal systems 

Kj = f (T,P) (4.13) 

For such cases, either direct curve fit expressions (Owens and 
Maddox^, 1968 ) or vapour pressure temperature correlations arc 
used. The most widely used expression for vapour pressure in 
the Antoine Equation, which, is of the fhrm 
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In (vapour pressure P ) =s a. 

VP j j 

where P^p j Is in mm of mercury and 
T Is in °C 


B 

1 * 


(4.14) 


Aj, Bj and are known as Antoine Equation constants and are 
characterisrtic of the particular compound under consideration. 
Extensive tabulation of these constants is available in Reid 
et al, (1977)* The partition coefficient K^. can be calculated 
from 


K. 


■ILJ 

p 


(4.15) 


^7here P is the system pressure in mm of mercury, 


4,3,2 Estimation of K^s Using Direct Methods ; 

The criteria that must be satisfied for equilibrium of 
a liquid and a vapour phase at temperature T and pressure P is 


fX = f J 1 = 1,2,'.,..C (4,16) 

3 3 

where f j is the f ugacity of component j in the vapour (V ) 
phase or liquid (L) j^ase and C is the number of components in 
the mixture. For a mixture where the vapour phase behaves 
like an ideal gas, one may write 



( 4 ’, 17) 


where yj is the mole fraction of component j in the vapour 
phase. The nonideality in the vapour phase, if any, can be 



taken into account by defining the fugacity coefficient 0j as 
follows 

Except for mixtures containing strongly associating 

V 

substances such as organic acids, the value of 0^ will be 
close to unity at low pressures. 

Similarly, the nonideality of the liquid phase can be ^ 

accounted for by considering a liquid phase fugacity coefficient^ 

Xj 

0j as follows : 

f j = 0j Xj P (4.19) 


where Xj is the mole fraction of component j in the liquid 
phase. 


Prom equation 4.16, at equilibrium 





X . 


J 


P 


y , j3 . . I 

or K. = — ^ (by definition) = -rl— (4,20) ! 

J 0j 

The above procedure for calculating partition coefficients 
has the limitation that a single equation of state is required I 
Which is valid for both the liquid and vapour phases for the 
whole rar^e of temperature and pressures. For this purpose, as i 
high as eight parameter equations of state have been proposed 
'’e.g;, the Benedict - Webb *- Rubin equation of state). Hbxvevcr, 
such procedures are highly uneconomic with respect to computatij 
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time and memory requirements specially for distillation 
calculations. Some of the three parameter equations of 
state (SRK or Peng*«Robinson) have been reported to give good 
results, particularly for hydrocarbon system and petroleum, 
fractions. 

For the calculation of partition coefficients by this 
method, the equation of state is used to calculate the 0^ and 

and the partition coefficient is calculated from equation 4,20, 
The Soave modification of the Redlich—Kwong equation (SRK)method 
and the Peng-Robinson equation have been presented in Table 3 
and 4 in different forms applicable to pure components and 
mixtures, . . 

It is convenient to rearrange Equation 4,46 by multi** 
plying and dividing the terms on the RHS by P/RT) and simplifying 



(4.21) 

(4.22) 

(4.23 ) 

(4.24) 



Equation 4.21 has three roots and the success of the 
single equation of state method is due to the fact that the 
highest real root can be taken to correspond to the vapour 
phase and the lowest real root to the liquid phase. The 
method for finding the cubic roots is given in Appendix E. 

Once Z has been calculated^ v can be found fran 

ZRT (4.25) 

V - -5- 


The fugacity coefficients required to calculate K^'s for 


both the liquid and the vapour phases can be evaluated by 
using the follcwlng derivation. If anS represent the 

ohemlcai potential and number of moles of the jth component, 

we know 




-P-i 






(4.26) 


(4.27) 


Again^ by definition for mixtures: 

= RT d In f j (at constant T) 

Integrating Eqn, 4.28/ 


(4.28) 


(4.29) 


Choosing the reference state to be ideal, we have 


— 

1 = ^3 


^4 P 
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Combining equations 4.27 and 4.29^ we have 


r , r c? 

J ^ 

£ . 

RT In —2- = RT In jZJ . 

tr R J 


(A-A 


•>] 

*«M#» rn 


(4.30) 


r,v, Uj. / j 


1-^. 


1 

Jt,v, 


(4.31) 


Again, if A is the Helmholtz free energy, we know that 

dA. = -PdV 


or AwA 


V 

- J 


P dV 


(4.32) 


Now using equation 4.42 and the SRK equatign of state- and then 
using the result in equation 4.31, finally, we have 


In ^ (Z-1) + in 2r^bT 


+ ( ^|j) In (“) *''-(TERM) 


(4.33) 


\;here 


(a .a,J’ 


TERM = ^ - 2 £ (1 - K ) 


(4.34) 


v;ith 


Zj^ = liquid phase mole fraction to give 0j\ and the 
vapour phase mole fraction to give 0^ 


RT bj Fj 


(4.35) 
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ilb 


RT bF 


(4,36) 


A similar procedure can be follov/ed to calculate fugacitp 
coefficients using Peng-Robinson equation of state in place of 
the SRK equation. 

The partition coefficient can thus be calculatSd by 
using equation 4,20. 

The results obtained for nearly ideal and nonideal systems 
have been shown in Chapter 5 (Problems 1 and 2). 

Using Indirect Me thods; 

In the indirect methods^ the liquid phase nonideality 
is characterized by the activity coefficient, defined as 


f ^ X . f° ; 

1 j J J 


j “ 1#’2,,,.»,C 


(4.37) 


V7here fj is the fugacity of the jth component at any arbitearily 
chosen standard state. Choosing the standard state as the 
pure component at saturation temperature and pressure, we have 


^o ofSat i-,sat -m/ -nSat 

f . = 0 . P . ... P . 

J J J J 


(4.38) 


vSat 


Since is close to unity at the saturation temperature 

of the pure component j. Therefore, from equations 4. IS, 4.18 
and 4.47, we have at equilibrium, 

/j Xj (4.39) 

i 




y. 


■*1 i 

or ^ = K, (by definition) = ^ 


(4.40) 


0\ P 
J 
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In these methods^ the term can be calculated^ as 
discussed earlier/ using an equation of state or correlations 
and the term Pj is obtained from the Antoine equation 
(Equation 4,14), 

It should be noted that equation 4,40 is structured in 

a way that the thermodyn anic functions have a minimam inter- 

V 

dejoendence on basic variables. Thus defends on vapour 

phase Composition/ P and T but not on the liquid phase 

conposition. Conversely/ ^ ^ is a function of liquid phase 

,sat 


composition only. The saturation pressure P. 
dependent on system temperature/ T, 


is only 




Of the various methods available to cc«ipute the 

Ul'JIFAC Group contribution technique is the most recent and 
widely used. In this method/the activity coefficients in 
the mixture are related to the interactions between the 
structural groups. The molecular activity cgefficient is 
divided into tv/o parts; a combinational part/ which is due to 
the difference in size and shape of the molecules in the 
mixture and a residual part which is essentially due to energy 
interactions amorig the functional groups in the mixture. For 
molecule j in the mixture 


In o . = In + In 

3 3 j 

Combina- Residual 
torial 


(4,41) 
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The theoretical aspects and detailed algorithm for 
calculation of and have been elaborately . dealt by 

I’redenslund (1977). Table 5 summarizes the final equations, 

4,4 Calculation of Enthalpies, h^^ . and Hu".: 

Another thermodynamic property, whj.ch is extensively 

used in distillation calculations, is enthalpy. For liquids, 

the enthalpy is abbreviated as h. . and for vapours it is H. 

i/j ■i/j 

For convenience sake, we will drop the subscript i in this 
chapter and denote them by hj and Hj, respectively. 

Curve fit experimental correlation for the enthalpies of 
all commonly used compounds are available in literature (Owens 
and Maddox 1968 etc,). Yen aand Alexander (1965) have given 
complicated generalised correlations for most of the compounds 
for various ranges of conditions. When estimation is required, 
the enthalpy departure function is calculated from the 
equations of state. The method ^f derivation for the departure 
functions is outlined below. These functions have been listed 
in Tables 3 and 4, It should be mentioned that they refer 
to mixture properties and do not have any subscript. 

For constant temperature, change in Helmholtz from 
energy A is given by 

m 

dA = -PdV (4,42) 

Integrating this at constant temperature from the reference 
state volume, V° to the system volume V, we have 
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^ V V _ % 

® P dV = - J" P<3V - J PdV 


V’ 


<30 

V 


V 


=-/ <P-P) 

<2Q 


V 


. d^^ - RT Itt 

V 


(4,43) 


Other thermodynamic functions can be similarly derived. 
The major departure functions of interest are 

V 


S-S" = « 


t)T 

V 


■ / MV-?] 


OCD 


dV + R In ^ (4.44) 

V° 


and 


= (A-A°) + T(S-S“) + RT (Z-1) 


(4,45) 


Thus, xvhen the equation of state to be used is selected,., 
the departure functions can be readily calculated* It must 
be noted that Eq,4,45 also applies to liquid enthalpies (h—h*^) 
’'.jfhen liquid phase mole fractions and compressibility factor 
are used. 


4,5 Calculation of Dew Point (T^^) and Bubble Point (Tg) 

The dew point temperature (T^), at any given pressure, is 
the temperature at which a vapor of given composition ,2 

gives .tiie first "drop" of liquid. Similarly, the bubble poi.nt 
temperature (Tg), at any given pressure, _ is the temperature at 
which a liquid of composition Z^, Z 2 , , ,Z^ gives the first 
‘bubble' of vapour. It has been mentioned earlier that dew 
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point and bubble point temperatures can be calculated from 
the follov</ing nonlinear equations by using a convergence 
technique (for example, Regula-»Falsi method), 

* z. ' 

fp = ^ -1 = 0 at dew point vrith ^j~yj (Qiven) 

j=i J 

%= 51 Kj Zj — 1 = 0 at bubble point xirith Zj=cXj (given) 

j=l 

However, it should be emphasized that in general, the Kj‘s 
depend on both the vapour and liquid compositions and therefore, 
the computation of dew points and bubble points require an 
iterative solution over two computational loops — one for the 
unknoxvn conpositions (vapour composition for bubble point 
estimation and liquid composition for dexi? point estimation) and 
the other for the required temperature. Tables 6 and 7 give 
the flol^r chart for these computations and problem 3 (Chapter 5) 
shows the tested example with the results. 



65 


TABLE 2 


Estimated Basic Physical Properties for Halo.iSubstituted Silane 

Compounds 


Pr opert ies 

D 

T 

M 

DH 

CL 

TH 

.■^e£ ^ . 


343.15 

339,15 

331,05 

314.15 329.95 

305,15 

ii 

M 

129.09 

149,59 

108.59 

115,09 170.09 

135.59 

- 

T^C’K) 

519.779' 

515.583 

499,659 

487.853 

503,436 

475.727 

e4. 

4.1 

P^(atm) 

33.128 

35.012 

30.691 

37,709 

36.456 

40.049 

Eq, ; 




. 



- 

4*? 

0\ G 

7.079 

7.111 

7,009 

6^866' 

7,115 

6,899 

Eci, 

4,5 

- 

344.189 

322.336 

358,901 

288,082 

302,193 

263.895 

Eq. 

"d(atm) (gmole) 






4.4 


0,2673 

0,26675 

0.26864 

0,2713 

0,26667 

0.2707 

e4 . 

4.6 


0.2657 

0.2722 

0.2513 

0.2219 

0.2729 

0,2287 

Eq, 

4.12 

i 


■'.;here D = Dimethyldichlorosilane 


T = Methyltriclolorosilane 
M = Trimethylchlorosilane 
DH = Methyldichlorosilane 
CL = Tetrachlorosilane 
TH = Trichlorosilane 

^Values taken iron private correspondence (Defence Scientific 
Information and Documentation Centre^ Delhi), 
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TABLE 3 


The Soave Redlich - Kwonq Ecruation of State 


For Pure Components: 

7 _ W V 
" RT t V-b ) 


Redlich-Kurong F = T^ • 


-H-a b „ 

ils 


(4,46) 


Soave 


F = i £i+( 0.48 + l,574(<?-0.176a>^) 


(1 - t; 


,0,5^J2 


0*4274802327; = 0,086640350 

For mixt ures ; (j refers to componant numloer) 


/-\ R Tc . 

= M- 1 . 

b Pc . ' 


>== Z 

d'l 


Tc 


j T 


P 


D 

-fi. 


480 + 1,574<>3j ~ 0,176 




J -fi-b 

c 

E 

j 


RT b . F . 
J J 


C C - Tc . P , Tc. F. -J 

T T \ * -K-] 


0.5 



Pc, 


lc.=l 


b . P 

Bj = -^“7 


B = 


e 


j=l 


Z . B . 

J J 
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a = jjij-S RTPb 7 

Jlh 


A = S pB 

jQ- b 


+ Z(A-.B-.B^) - AB :i= 0 


(4.47) 


inp. - ^ (Z-1) + In 


V , a T^/V+b^* 

^-b7 bRT v “^ 


i 


(l-K^-j^) 


(a. a. ) 


(4.48) 


Zj = liquid phase mole fraction to get 0j 

. J7 

= vapour phase mole fraction to gat 0j 
Expressions for the Departure Functions; 


energy A: 


A-A'’ = RT : -ln(^) - ^ ln(— ) - In(^) 

o 


v;here 


(4.49) 


Y* "S” Z,Z„ (l^K ..)/”* (^^ )*(^^)b .b,F F 


j-1 K=1 


(RT/P° ) 


Molai 


= R (In (^) - 1“ * In(^) + in ^^) (4.50) 


^ 21 2 [j- 


v/here 


b .b, F .F, 
J k 3 3^ 


j=l k=l 


r-f(<^j^) f(<^.) 

I "To^ 

^ k J 



68 


f(^j) = 0,480 + 1,574 - 0,176 



H-H° = (JWA") + T(a-S®) + RT(Z*1) 

T - • 

= J" Cp dT (i,e. Zero at standard condition 
^o chosen, T°) 

P" = 1 atmosphere. 
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TABLE 4 


The Peng - Robinson Ecfuation of State 


For Pure Compo.nents ; 


RT 


(v— b ) v(v+b) + b(v— b) 

where a and b are constants 


(4,51) 


For mixt ures : (j denotes component numbers) 

RT 

b . = 0,7780 




b = y 2 b . 

J J 


C,J 


j=l 

0,47724 


r,2 


0,5 


cCj, = 1 + (0,37464 + 1.54226 6A-0*26992 65p(l-T°'^) 
J J J j 


C G 

j=l k=l 


0.5 0,5 


where K., *s are the binary interaction parameters (They are 
taken to be unity where sufficient data are not available). 


B =5? ' 


. a P 
A = -2 — 2- 
R 

2 *_^„ ^2 


2^ +(B«1)2^ + (A-.3B^-2B)2 - (AB-B -B^ ) == 0 (4,52) 
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Table 4 (contd) 


V * ZRT/P 

ajj "" Sj * 

b 

In JZfjw ^ (Z-.1) - In(Z-b) 


2;329B 


— 


■J 


I In f 2^'2*414B x 

a “ b / ^z-^),414B 


(4,53) 


•where Zj » liquid composition to get 0j and 

V 

= vapour composition to get Pj 
Expressions for Departure Functions for Mixtures : 

.0 


£( o5.) = 0.37464 + 1.54226 o9. - 0.26992 OS 


c . c 


da 

dT 


' - ir £ X (0.45724°-5r) 

i=l k=l I- 


J 

tO. 5 f ( 
^k 




,0,5 




tO.5 £ ( ^3 ) 

!iL_ L ^ 

„0.5 k 


da 


-j 


V - 


s*-s^ 


O = ( )in 2414 R In (— 2>^) (4,54) 

^ 2 b’ V. «. n 414 b' V — b 


Z - 0.414 B 


V — b' 


where, S denotes entropy of the mixture 


da 


(Z + 2,414 B) 

H-H° = RT(Zwl) + j' • In 

K2 si ^ (Z « 0.414 B) 


(4.55) 


vjhere, H denotes enthalpy of the mixture* 
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and A-,A° * - T(S-S®) - RT(Z-l) 

C 

dT 

T ^ 

298°K 

T 

H° » J Cp cfT 
298 °K 

v° = RT/P° ; = 1 atmosphere 
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TABLE 5 


Uhifac Group Contribution Eq-uations 


In 


= ln-" = 


3 ^ 3 

Combinatorial 


In 


>J R 


(4.56) 


Residual 

0 . 


In ^ + I ^ S' ^ ^ ^ 3^3 


j=l 


xvhere 


Ij =-|(rj-qj) - (rj-l); 2=10 


q . X . 

nil- , 


S ^ c 

j=i 

Molecular surface 
area fraction 

The vender Waals volume r 


0 


r . X . 
^ L 


3 c 

“Z 

j=l 

Molecular volume 
fraction 

< (i) 

k=l 

N 

(i) 


The vander Waals surface area q^ = r Si: 


Q, 


k=l 


where N = number of groups in molecule i, 

Rj^ and Qj^ are group volume and area parameters tabulated in 
Fredenslund (1977) 


In 




R 


k 


5"(i) 


k in r k - 1 " I'' k 


(i) 


(4.58) 


all groups 

K = 1,2,...,N (where N = nurriber of different groups in the 


mixture ) 
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\ 



= residual activity coefficient of group k in a solution 

= residual activity coefficient of group k in a reference 
solution containing only molecules of type i 


In r\. = Qj^ , l-ln/X" j^mk^ “ ^ ^ km/ ^nm^ (4.59) 

^ ^ n 


m andn = 1,2^.*..,N (all groups) 
Eqn*4*54 also holds for 


(i) 




(P) 


o = 

m 


Q X 
m m 

N' 

Q X 

/ , n n 
n=l 


Group surface area 
fraction 


m 


X 


X 


- 


m 


p=l n=l 
Group fraction 




= exp (*-a ,/T) 

nm nm 


(4,60) 


where is the binary group interaction parameter given in 

Fredenslund (1977). It should be noted that there are two such 
parameters for each pair of groups and a^^ / a^^ 
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TABLE 6 

Pldwchart for Dew Point Temperature (T^) Estimation 


START 


INPUT DATA ; C ^ 2 . ^ ; EPSLON , 
P,S(Usually^0,6 ), 
initial guess for T^ 


j “ j' ""J/Old 
I COUNT = O 
K FLAG = 1 





^J,neiv'"'^j,old 1 
"^ss than EP^ 
\(for ai;L^^ 




• j , old^ j7new ; ICOUNT=ICOUWT+ 


Ca leu late new esti- 
mate for dew point 
Tj^ using Regula falsi 
technique 



PRINT 

I?„ 


■ 




STOP 
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TABLE 7 


Plow Chart for Bubble Point Temperature (T^ ) 

— ^ „ x5 

Estimation 







CHAPTER 5 


RESULTS AND DISCUSSION 

i^veral test problems have been solved using the 
algorithms developed during the present work. Through these 
examples an attempt has been made to check the programs for 
VLE data prediction, basic thermodynamic property estimation, 
dew point and bubble point estimations, distillation column 
simulation and design having numerical schemes like Newton- 
Raphson, Regula-Falsi, Thomas Algorithm, Block ’^omas Algorithm, 
Gaussian Elimination, Golden Section Search etc. Some of 
these problem categories will be presented in this chapter. 

They are as follows ; 

Problem 1; Estimation of K~values for nearly ideal and 

I 

nonideal systems 

Problem 2: Binary VLE data prediction for Halogen-substituted | 
silane systems j 

Problem 3 : Dew point and Bubble point calculations for ideal 
and nonideal systems 

Problem 4: Distillation Column Simulation using ideal and 
nonideal thermodynamics 

Problem 5: Distillation column design using the new simultaneous 

-f - - , [, 

semi— tray-by— tray formulation, 

[: 

Comparison has been made of computed results xi/ith published data I 
wheravear possible which validates the developed programs. ; 
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VLE Data Prediction (Problem 1 and 2) : 

Problem 1 demonstrates the estimation of K»*values for 
common systems for which all the critical values are available 
Both Soave-Redlich-KWong (SRK) and Peng— Robinson equations 
of state have been used for prediction of the K-values, 
Available experimental values which are used for comparison are 
quite reliable. The agreement between the two is generally 
good. The follaving observations can however be made: 

(i) the accuracy for prediction of VLE data by 
equations of state decreases at conditions 
near and above the critical point (e.g, 
methane in problem 1,1), 

(ii) the accuracy of prediction also decreases 
for Compounds having associations (e.g. 
methanol in problem 1.3), This observation 
has also been made by many other authors _ 

(Predenslund^ 1977, Srinivasan, 1982 etc.) 

Problem 2 presents the prediction of K— values for systems 
(Halos ubsti tut ed silanes) for which critical properties are 
not available and have been estimated by the procedure 
outlined earlier. The input data are molecular weights and 
boiling points of the various compounds. Only meagre experi- 
mental data is available on these systems and that too vjith 
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limited, reliability. The agreement betv/een predicted results 
and experimentally observed values is gen car ally good except 
in case o£ silicon tetrachloride and trimethylchlorosilane 
(at certain conditions) for which the equation of state may be 
less applicable or the observed values may be less certain. 
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Ijj; Estimation of Partition Coefficients (K-values) for 


Nearly Ideal System 
System: A mixture of Methane 

Ethane 

Propane 

Temperature : 283,15 K ; 

= 0,029 ; 

X2 = 0,166 ; 

X3 = 0,805 } 


( 1 ) 

( 2 ) 

(3) 

Pressure: 13,61 atmosphere 
= 0,290 
= 0,292 
Y2 = 0,418 


Physical Properties 


Compos 

nents 

M 

*^0 

(■’K) 

P 

c 

(atm) 

V 

c 

( cm^/g .mole ) 

Z 

c 


1 

16,043 

190,6 

45,4 

99,0 

0,288 

0,008 

2 

30,07 

305,4 

48,2 

148,0 

0,285 

0,098 

3 

44.097 

369.8 

41^9 

203.0 

0,281 

0.152 


Comparison of Kj values : 


Component 

Predicted values using 

Experimental 

values* 

SRK Eq, 

HHESHHHHii 

1 

9,2 

9,0 

10,0 

2 

1,82 

1,8 

1,76 ; . ' 

3 

0^545 

0.543 

0.52 


Taken from Smith and Vanness (1975) 
CPU time: 0.07 sec 
Elapsed time: 0,76 sec. 
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Problem 1,2: Estimation 

of 

Partition Coefficients for Nonideal 

System: 



System: A mixture 

of 

Benzene (1) 



Cyclohexane ( 2 ) 

. 


n-hexane (3) 

Temperature: 342,7 K 

« 

/ 

Pressure: 1 atmosphere 

= 0,279 


y^^ = 0,250 

X2 = 0,037 

/ 

Y2 - 0,025 

Xg = 0.684 

; 

y^ = 0,725 


Physical Propectles ; 


Canpo«» 

nent 

M 

T 

c 

(°K) 

P 

c 

(atm) 

Vc, . 

3 

(cm /g.mole) 

2 

C 


1 

78,114 

562.1 

48,3 

259,0 

0,271 

0,212 

2 

84,162 

553,4 

40,2 

308,0 

0.273 

0,213 

3 

86.178 

507,4 

29.3 

370.0 

0.260 

0,296 


Comparison of K. Values: 

<m - p,,.. - - 


Component 

Predicted 

Cj values using 

Experimental 


SRK 

Peng .Rob in son Eq 

values* 

1 

0,8445 

0,8605 

0,8960 

2 

0,7442 i 

0,7551 i 

0,6756 

3 

1.0375 

1.0397 

1,0599 


Taken from Srinivasan (1982) 
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Problem 1,3 : Estimation of Partition Coefficienbs for a Highly 
Nonideal System ; 

System: A mixture of Methanol (1) 

Carbontetrachlor ide ( 2 ) 

Benzene (3) 

0 

Temperature: 328,15 K ? Pressure = 0,9296 atmospheres 


^1 

= 0,194 

7 

^1 

== 0.507 

X 2 

= 0.592 

• 

^2 

= 0.373 

X 3 

= 0.213 

• 

/ 

’^3 

= 0.12 


Physical Properties 


Compo- 

nent 

M 

T 

c 

(“K) 

P 

C 

(atm ) 

V 

3 C 

(cm /g.mole) 

2 c 




1 

32.042 

512.6 

79,9 

118.0 

0.224 

0.559 

2 

153.823 

556.4 

45,0 

276,0 

0.272 

0.194 

3 

1 

78.114 

562,1 

48,3 

259.0 

0,271 

0.212 


Comparison of K, values; 


Component 

Predicted K 

j values using 

Experimental 

Values*^ 


SRK Eq, 

Peng “Robins on 

Eq. 

1 

0.97 

1.015 

2.613 

2 

0,5462 

. 1 

0.5601 

0.63 

3 

0.4697 

0.483 

0,5634 


Taken from Predenslund et al, (1977) and Srinivasan (1982) 
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Problem 2,1 ; Binary VLE Data Prediction for Halosubstit uted 
Silane Systems 

(a) System: A mixture of Trlmethylchlorosilane (1) 

and Methyltrichlorosilane (2) 
Pressure: 1 atmosphere 

Physical properties: taken frcm Table 2 
Conparison of K. values: 


Conditions I 

Compo*- 

nent 

Predicted K . values 
using-^ 

Experimental 

values* 

Tempera 

Mole 

ture ( "ic) 

fraction 

SRK Eq. 

Peng “'Robin- 
son _ _ _ 


x-so*i 


1*2446 

1*2455 

1*19 

338,05 

yi=olyll9 

2 

0,9443 

0,9499 

0.9789 

337,25 

x^=0i25 

1 

1.2046 

1*2059 

1*18 


Y 2. ^ ^5 

2 

0,9229 

0,9288 

0.94 

335.35 

x^=0,50 

: 1 

i 

1 

1,1266 

1.1288 

1,136 


y2=0;568 

2 

0,8754 

0.8817 

0.864 

333.05 

Xj^=0475 

1 

1.0454 

1*0486 

1*0747 


y^^ =0,806 

2 

0.8221 

0.8290 

0,776 

331.57 

X^=0*9 

1 

1 

0,9980 

1*0017 

1,0288 


y^^ =0,926 

2 ■ 

0.7900 

0,7972 

0,74 


*Taken from Zel'vensky (1971) 










(b ) System: A mixture of Trimethylchlorosilane (1) 

and Dimethyldichlorosilane (2) 
Pressure; 1 atmosphere 
Physical properties ; taken from Table 2 
Comparison of values 


Conditions 

Compo- 

nents 

Predicted K . values 
using ^ ' 

Experimental 

Values* 

Temper- 

Mole 

ature 

(‘‘K) 

fraction 

SRK 

'•Eq 

Peng- 

Robinson Eq 



“~I 


liOlsT^ ' 

1.33 

342,15 

yj^ =0,135 

2 


0,9549 

0,9611 

340,35 

x^=0*25 

1 

1,0145 

1,02 

1,248 


yi =0,312 

2 

0,8993 

0,905 

0,9173 

337.17 

x^=0,5 

1 

0,9189 

0.925 

1,2 


yi=o,6 

2 

0,815 

0,8213 

0.8 

333.93 

Xi=0,75 

1 

0.829 

0,8356 

1,1026 


y3_=0,827 

2 

0.7357 

0,'74 ■ 

0,692 

331.95 

x^isi0,9 

1 

0,7775 

i 

0,7846 

1,0411 


y^^ =0.937 

2 

0,69 ; 

0.697 

0,63 


if 

Taken from Zel'vensky (1971) 
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(c) System: A mixture of Methyltrichlorosilane (1) 

and Tetrachlorosilane (2) 

Pressure: 1 atmosphere 
Physical properties: taken from Table 2. 
Comparison of Kj values 


Conditions 

Compo- 

nents 

Predicted K. values 
using ' 

^ 

Experimental 



values * 

pRBSH 


mipggiiiii 

ESS 

- 

337,05 

, 1 ,-" 

X^=0-,9 

1 

0.9149 

- ' — -Z.Z — ' 

0,9209 

0,95 


y]_=0.855 

2 

1,2151 

1.2184 

1.45 

335.09 


1 

0,8601 

0,8665 

0,9013 


y2=0}676 

2 

1.1448 

1^1489 

1,296 

332.62 

x^=0.5 

1 

0,7946 

0.8015 

0,866 


y ^4 3 3 

2 

1 ^0606 

1.0655 

1,134 

331.30 

x^=0,25 

1 

0.7612 

0.7684 

0,844 


y3_=o'“211 

1 2 

1,0175 

1.0229 

1,052 

- 

X. =0,1 

1 

0.739 

0,7464 

0,83 

330,40 

1 X 






y ^083 

i 

2 

0,9889 

i 

0,9946 , 

1,0189 


ic 

Taken from Zel'vensky (1971) 










5*2 Dew Point and Bubble Point Calculations (Problem 3 ) : 


Calculations have been made using ideal as well as 
nonideal thermodynamics. In problem 3,1 (a) and 3,2 (a) ideal 
thermodynamics is used to estimate the dew point and bubble 
point temperatures of vapour and liquid ^^^hich are themselves 
in the thermodynamic equilibrium with each other. As shewn 
below, the results are in close agreanent with publi^ed 
values (Amu- . , 1958). In problems 3.1(a) and 3.2(b) the 

system is same as that inpart (a) except that nonideal 
thermodynamics is used to calculate K-values, This, however, 
changes the mole fractions of various components in the two 
phases which are in thermodynamic equilibrium. The comparison 
of results with part (a) shows hew these temperatures vary 
when nonideality is introduced. Problems 3.1(c) and 3.2(e) 
use the same vapour and liquid compositions as in 3. 1(a) and 
3,2 (a) and proceed to calculate the two temperatures using 
nonideal thermodynamics. These compositions are different 
fron these in part (b) where the compositions of vapour 
and liquid themselves were calculated using nonideal thermo- 
dynamics. 

Problem 3.1(a) 

Feed composition: 

n-C^d) s 139,4210 lb,moles/hr 

n-C^(2) : 12.7402 Ibjmoles/hr 

I 0,7588 X. I0‘"^lb.molos/hr 


n— Cg ( 3 ) 
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Thermodynamics used : Ideal 
Initial guesses supplied: 120'’P and 170*P 
Calculated dew point temperature: 144,044'*F 
Literature values : 140®P 

CPU time : 0,07 seconds 

Elapsed time : 3,10 seconds 

Tolerance limit : 0,0001 

Problem 3,1 (b) 

Feed composition: 

*4 * 

n-C^ (1) : 121,739 lb,moles/hr 

n-C^ (2) : 14.1454 ,, 

n-Cg (5) : 0,20218 ,, 

Thermodynamics used : Nonideal 

Initial guesses supplied; 120®F and 170‘’P 

The liquid mole fractions 

estimated : (1) 0,773985 

(2) X2» 0.2195954 - 

(3) X 3 = 6.41947x10*^ 

Calculated dew point temperature; 140,64*P 

Number of iterations : 4 

CPU time : 1,47 seconds 

Elapsed time ; 27,56 seconds 

Tolerance limit ;; 0,0001 

Problem 3.1(c) ’ 

Feed composition : Same as in 3,1 (a) 

Thermodynamics used : Nonideal 

Initial guesses ; same as in 3,1 (a) 

No, of iterations s 4 



Calculated dew point 
tfaonperature 


: 137.87«P 


CPU time 
Elapsed time 
Tolerance limit 
Problem 3.2(a) 

Feed composition: 
n-C^d) 
n-C^(2) 
n-C5{3) 

Thermodynamics used 

Initial guesses supplied 

Calculated bubble point 
temperat ure 

Literature value 


; 1,46 seconds 
: 29.82 seconds 
: 0.0001 


: 63,7496 lb .mole s/hr 
: 40.077 ,, 

: 2,5466 ,, 

: Ideal 

: 150®P and 180®P 

: 163.7533‘’P 
: 160. 53 


CPU time 

Elapsed time 

Tolerance limit 
Problem 3.2(b) 

Feed Composition: 


n-^3(l) 

n~C^(2) 

n-C5(3) 



Thermodynamics used 
Initial guesses supplied 


: 0,05 seconds 
: 2,98 seconds 
: 0,0001 

: 46,7383 lb,moles/hr 
: 45.6364 lb,moles/hr 
: 4,9448 lb, moles/hr 
: Nonideal 
~ 150®P and 190®P 


The vapour mole fractions 
estimated 


Yl = 0,6871599 


Calculated bubble point 
tenperature 


72 = 0.2964676 

73 = 0.0163724 


170. 6868 °P 



Number of Iterations 

CPU Time 

Elapsed time 

Tolerance limit 
Problem 3, 2(c) 

Peed Ccmpos it ion 

^hemodynamics used 

Initial guesses 

Number of Iterations 

Calculated bubble point 
temperature v 

CPU time 

Elapsed t ime 

Tolerance limit 


: 3 

: 1^02 seconds 
: 20,56 seconds 
: 0.0001 

• same as 3, 2 (a) 

• Nonideal 

K 150®? and 180®F 
: 2 

! 157.3949®F 

• 0,62 seconds 
5 2.52 seconds 
5 0.0001 


5,3 Simulation of Distillation Colimns (Problem 4) ; 

The multiconponent ^ multistage distillation column 
calculation package has been used in simulation (or rating) 
mode without any multitray models to check its performance. 
Problem 4,1 deals with a 3 component feed to a 16 stage column 
(including reboiler and condenser). Both vapours and liquid 
are assumed to behave as ideal fluids. The calculated 
Component flowrates of distillate and bottor^ match very well 
with the publi^ed results (itmiundson, 1958), 

Problem 4,2 uses the same feed as in Problem 4.1 and the 
same, column conf iguration^but both vapours and liquid are 
treated as nonideal and the single equation of state appro^fh 





is used to account for the nonideality. The component flow 
rates obtained do not vary significantly fran those in Problem 
4,1 primarily because the feed mixture is near ideal. However^ 
the temperature prof iles in the two problems vary to a larger 
extent than the compositions (see Appendix 6), 

Problem 4,3 has been used with same feed conposition as 
before but. with a different column configuration (13 stages) 
Ideal thermodynamics has been used. As csxpected/ the separation 
achieved has decreased in this case as compared to problem 4,1, 


Feed rate 

Feed composition 
(component flew rate) 

Feed condition 
Feed temperature 
Reflux ratio 

Condenser temperature 
Condenser type 
Probl^ 4,1 Number of stage 
Feed location 


1CX),0 Ib,moles/hr with 

(1) n-C^ * 23.0 lb. moles/hr 

(2) n—C^ : 37,0 

(3) n^Cg : 40.0 
Saturated liquid 
225, O^F 

5 

22,6 lb, moles/hr 
120,0‘’P 

Total condenser 
16 

8th tray from the top 


Ibe following top and bottom flow rates were obtained using 
ideal thermodynamics. 
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Component Distillate, v, 

- - - - - .... j 


Bott ons, l^g 


21,5360 lb,tnoles/hr 1,4637 Ib.moles/hr 


n-^4 


1.0614 


35.9385 


0.0022 


39.9978 


nTotal 


22.5996 Ib.moles/hr 77,4000 lb .moles/hr 


Optimal values of ^ (in modified Newt on-^aph son method) were 


found to be between ,98 and .9999 for different iterations, 


Tolerance limit 


: 0,0001 


No, of Iteration 


CPU time 


: 4,28 seconds 


Elapsed time 


s 30,58 seconds 


The vapour and liquid flow and temperature profiles obtained ini 


the:' column have been provided in Appendix G, 


Problem 4,2: Number of stages : 16 


Feed location 


: 8th tray fron the top 


The following top and bottcm flow rates were obtained using 


nonideal thermodynamics. 


Component 


Distillate, v^^ 


Bott cms 1 


16. i 


21.47325 Ib.moles/hr 1,5267 lb,moles/hr 


1.11868 


35.8813 


n^C^ 


0.00806 


39,9919 


Total 


22,59999 Ib.moles/hr 77,3999 Ib.moles/hr 
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qptimal values of p (in modified Newton-Raphson method) 
were found to be between 0,94 and 0,999 for different 
iterations. 


Tolerance limit : 0,0001 

Number of iterations : 5 

CPU time : 25,54 seconds 

Elapsed time : 1 minute 29,52 seconds 

The output listing for the vari^le profiles obtained in the 
column has been given in Appendix G, 

Problem 4,3 : Number of stages: 13 

Feed location : 7th tray from the top 
The following top and bottom flow bates w(^e obtained using 
ideal thermodynamics. 


Conponent 

n-C3 

n-Cg 


Distillate, Vj^ 


Bottoms li 


20,8177 lb. moles/hr 2,1888 lb,moles/hr; 
1,7796 lb,rooles/hr 35,2203 lb, moles/hr 
0,0092 lb,moleSy/hr 39,9908 Ib.moles/hr 


Total 22,6065 Ib.moles/hr 77,3999 lb,moles/hrj 

W IIII 11 - ■ I.. » — ■■ ■■ m mj 

Optimal values of ^ (in modified Newton-Raphson method) were 
found to be between 0,98 and 0,9999 for different iterations. 
Tolerance limit : 0,0001 

No, of iterations s 3 I 

CPU time : 3.96 seconds 


Elapsed time 


: 35,8 seconds 



The output listing for the variable profiles obtained in the 
column has been given ini Aj^endlx Gf; 

5,4 Design of Distillation Columns (Problen 5) ; 

The earlier problems 4,1 and 4,3 have been solved in the 
design mode l.e,, where the number of stages and feed locations 
have been treated as unknowns. The vapor flow rate of component 
1 on unit number . and vapour flew rate of conponent 3 on unit 
number were used as design specification. In problem 5,1, the 
total number of stages was found to be 16,04 to be compared 
with 16 trays (problem 4,1) and feed location at 8,2 to be 
compared with 8, The Ednister approximation seems to be 
quite satisfactory for this problem as evidenced by a very 
small variation obtained, as compared to actual one. 

In problem 5.2 the total stages obtained is 13,02 to be 
compared with actual number of 13 plates and feed location 
at 7th plate exactly same as actual. 

Total feed rate : 100,0 Ib.moles/hr with j 

component feed rate (1) : 23^0 lb, moles/hr s 

(2) n-C^ : 37,0 lb,moles/hr 

(3) n-Cj- : 40,0 Ib.moles/hr | 

Feed condition: Saturated liquid I 

Peed tenpperature : 225, 0®F | 

Thermodynamics used : Ideal 

Sp€5Cif icat ions : 

Reflux ratio : 5 ; 

Distillate rate : 22,6 Ib.moles/hr I 

Condenser temperature: 120"P I 
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Pr pblgiiTt 5 y 1 i 

Vapour flow rate for component 

1 on unit 3 : 130.42100 Ib.molesAir 

Vapour flow rate for component 

3 on unit 9 : 17.374790 IbJjnoleSy/hr 

“rhe following top and bottom flow rates were obtained. 
Component Distillate^ ^ Bottoms^ ^11^1 

n -^3 21.5383 Ib.molesAir 1.4617 lb .mole s/hr 

n-C^ 1,0590 Ib.moles/hr 35,9409 lb .moles/hr 

n-u:g 0.0027 lb,moles/hr 39,9973 lb, moles/hr 


Total 


22,6000 lb, moles/hr 77,39999 lb,moles/hr 


Total number of trays • 16,04 

Feed location at 8,2 th tray from the top; optimal ^ was found 
to be between 0,96 and 0.9999 

Tolerance limit : O.OOOl 

Number of iterations : 5 


CPU time : 20,62 seconds 

Elapsed time : 58.5 Seconds 

Problem 5,2 : 

Vapour flowrate for component 1 on - 

unit 3 : 123,03 lb,moles/hr 

Vapour flow rate for component 3 

on unit 9 : 18,8633 Ib.moles/hr 

The following tojo and bottom flow rates were obtained: 



Component 


n-C- 


Distillate.v 


l,i Bottoms/ 


n-Cc 


Total 


20,8114 lb, moles/hr 2.1086 lb, moles/hr 
1,7794 lb, moles/hr 35,2206 lb, moles/hr 
0,0092 lb, moles/hr 39*9907 lb,moles/hr 


22*6000 Ib.moles/hr 77,3999 lb, moles/hr 


Total number of trays : 13,02 
Peed location at 7th tray from the top.. Optimal 
found to be between 0*96 and 0,9999 


was 


Tolerance limit 
Number of iterations 
CPU time 
Elapsed time 


: O.OOOl 
: 5 

: -20,62 seconds 
: 50,5 seconds 



CHAPTER 6 


CONCLUSION AND RECOMT^ENDATIONS 

In the present work, multistage, multicomponent 
distillation column design and simulation procedure has been 
developed* The existing literature lacked a gerffiralized 
doubly acting algorithm which could be applied to various 
column configurations and take into account nonidealities with, 
respect to thermodynamics and plate efficiencies*. The 
proposed formulation has been found to overcome most of the 
shortcomings* A combination of tray— by-tray and shortcut 
methods used in the present study enables the user to design 
a column of any desired configuration without repeated 
application of simulation procedure as previously required* 

The shortcut method applied to the multitray models has been 
made as rigorous as practical and most of the simplifying 
assumptions have been removed* The nunerical aspects of the 
formulation have been selected and Implemented in away that 
computation time and memory requirements are minimum. The 
computer package has been written in FORTRAN language and 
implemented on DEC— 1090 system* Some problems from literature 
have }Deen tested and good convergence characteristics have been 
obtained. The CPU time and memory requirements are almost 
independent of the number of trays. 
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It was found that about 60 percent of the CPU time 
required for the computations was taken up by the numerical 
Jacobian calculations. Instead of evaluating the Jacobian 
numerically in every iteration, and alternate procedure 
(Kaufmann, 1983 ) of updating the Jacobian from the previous 
iterations can be experimented upon. However, it should be 
recognized that using the previous iteration Jacobians will 
increase memory requirement since in that case it will have to 
be stared in the monory which in the j>^sent method is not 
necessary. Extensive variable packing schemes could be used 
to reduce memory requirements whicii have not been attempted 
in the existing program for the sake of clarity, of understanding. 
The effective stripping factor of Edmister may be 
treated as a funct:x.on of the number of stages inside the 
multitray unit, when its size becomes too large (Ohmura, 1978 ) 
or the size of the multitray unit may be shrunk with a corres- 
ponding increase in the single tray unit. The thermodynamics s 
for associated ccrapounds and highly nonideal systems and j 

property estimation near critical conditions should be further | 

1 

studied and improved upon. Calculations for e variety of I 

operating industrial distillation problems are required to be 
tested and checked against some reputed simulator packages (e,g,, ; 
Chemshare) regarding performance, memory and CPU time. This j 

should also build confidence in the present package with respect 

' " I 

to its design capability* ! 
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APPENDIX A 

Chart for Lev/ls~>Matheson Trav»*bv*.Trav Alaorithm (1932 


Specifications: feed informations, R, split of Key 
components, type of condenser and reboiler 

Iterative Variables: D, 


Assume d-, and i=l 


i/J 


Repeat a similar calculation 
starting from the bottom and 
check if (criteria 1) is 
satisfied. 


If (Criteria 1) is satisfied 
calculate corrected 

' / j 

by using Bonner's Eqn, or 

some other criteria 
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Specifications: Peed information, key conponent split, R,D,P, 

type of condenser and reboiler 
Assume: T , 1. ., v. ., N„ and Ne 



Inner loop of Najiitali^Sandholm iterative 
solution to obtain corrected T, , 1. - . and 

V. . ^ 

1/1 


Estimate changes in (N_+Ng) and (Ng/Ng ) 


/change\^ 

in Ng and 
No 0.5?^ 


Ad just 
Ng and Ng 


Estimate new T., 
: V. . and 1 

j j-/ j 

' by scaling 


'all f unct- ^ 
ions 4=€ ^ 


Print the results and 
output variables 
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APPENDCX C 

Flowchart for Ohmura Kasahara Algor itlim (1978) 


START 
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APPENDIX D 

M Chart for the New Simultaneous Semi—tray—b'' 


START 


Read input variables; 
specify parameters 



Use Block Thomas Algo- 
rithm to get ^correct ion 
vector, 


Calculate initial guess values 
for and T- 



Calculate component liquid 
and vapour flowrate from 
tridiagonal matrix solution 
of combined material balance 
and equilibrium relations 


Use Golden Section 
Search to get Optimal 


Check 




olersncc 


limit 


Calculate discrepancy 
functions: P 


Print output 
and result 


Calculate Jacobians ; J 
(numerically, analytically 
or mixed mode) 


Evaluate new v-'ri* 
able ve'ctor ^ ’ 

new old / 
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APPENDIX E 


Method to Obtain Solution of Cubic Equat ion (Perry, 1963 ) 
3 2 

Equation: x' + tax + cx + d = O 
p = (3c ^ b^) 

q = (27d^ - 9bc + 2b^) 

R = (p/3)^ + (q/2)^ 

M = (-q/2 + /5 positive value of R used 
N = (— q[/2 — dO“ 


Roots 

for 

r'>^o 


= M + N 

72 = (M+N) + 

^3 = - | _(M+N) ^ 



^3 


2 

i 


(M-N) 

(M-N) 


i = 1,2,3 


If 

If 

If 


Roots 
for 
R Co 


R one real two cofnplex conjugate roots 

R = 0, three real roots, of which, at least tv7o 
are equal 

R ^O, above formula' are impractical, then use 

= (+ 2 y^-p73 ) cos PI + 120°K*| - K = 0,1,2 
2 

0 = Cos"*^ ( ■ : — ~), positive root taken 
-p /27 

Use upper (-) sign for qj>0, Icwer (+) sign if q^CP 
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NsNV 


PACK) .K?:j , J)/-l ,OUV. 

PiiCK)?Ksl?3)/7 Jl5E-02.7.f t3G-02.Kl|4|-ul/, 

.PC ? K=1 . 3 ) /- 3 , 789E-05 . -2 ,b--i 7b-0b , -6 , 5 Uh 

.Pfl. cK } 1 3.) /.I ,, 91 8E--09 r -0 • l> 1 4E-«i39 ^ 1 • 2 ^ ^ 

‘■i^ODefiOLSfi data fur KOfUDeAU THfcRKODY:'iAATvS 

(PCCK) ,K.=i .3)/41i9.37.5,33.3/ . , 


5/ 


Uv.aL 7 .TEu 

■ \f * -w 4 . 4 ^' - 4 . *,.' » ' '••-* x. *-* *N( ' # S. »* *.»' ; 1 X T * Tl Tv * I* . 

UtOUID Er4THAliPY IS OHuY CkhOHAl'^V 


InOEX= 1 ? ^EbUILIBBIUN CALCiJLATTQN IS DD^> 
I? 40 eXs .3 ? LtlQ 


0»f SlPT7 / 1 4 « 7 

Ts(CT|-12lO)n5.0/9.O)) + 27 3.l5 
IFCIHDEX.GE-2) go to 275 

tRasp-eoa? heturh 

SOfU^O.O 

SUi2»0,0 

00 126 X=lfN, ^ 
SJ|«l=SUMH-XX(JvSP,K) 


iimr 


) 


JUKI 

DO 127 K=N+1.2^H 
S0M2sSnK2+XXCaSP.K) 

COIJTTMUS 
YAP = .5U'12 

DO 12B K=t .0 , , 

X(F)=XXCJSP.K)/RfaO 
Y(K3 = XX(.TOP,K+^OXVAP 
cij'ii'iuaE 

C A U U c. U i ! L i<i ( t » P f fXfY/KI) 

DO 

RKCJAt'.’O'sElCK) . 

cn-M'T^uv. 

IF (.ip-^T TTU.D) GO TO SOT 
If ,LK;>XTO,e:v-! 3 GO TO 20 

F[}iVSrri4Vi^EULa'* corat by pfkg-robxmso,'^ eo^.') 
r ,0 };u . )*' 

FOP V^TCT^i ^ CvOjAM COOST by snAVF>REOIjICH*KWuNG EQW 
CO'Mi'li'iUe: 

f"! f 3 4 4 1 V'l ^ 1 - # a 

pSL«TcTi-eT’3l’-=04UJxO»'If'*) = -,E16.8,4X,-rO|! COMP. Nu. 

cu-n'iro.!-; 

r.n It) 

XFn.'iuEX*eD. )) GQ.TO 27ti 
IS'hlSO, T). i ) 

DO IflrK. 

DO 27o !jL = 1 f < 

YDO'Uua-o.*-- ... . 

If Ci.bO.UO) tl5'J-H i)r1 ,0 

CAUpaioTf ZnAX,XMIt,yi)UM,F,B.FM,TC, PC, T.PT, OMEGAS, H.KSITU) 

PCS)AT = i>CCi) 

TCl>Ar='rca) 

2CDAT=ZCC I) , 

CPAxcpv.tPrtti) : 

cpBscpy.Hptri ) / : 

. CPC=C?V\?CCI ) 

81m/ MT oil ( Pt , r , Z ,?CDAT, PCDAT .ZCDAT ,CPA ,CPB ,CP 

i CP^HAT) 
cpFCjjjpa jfCPvoM 
fe>)TYCJv>i’.U = a;*CC 
DEI.HYCDsOFiVt , 


M3) 


■•'C,GP 0 ,riACT,DFl,t 4 


DFSJI^/: 


i'H; 

GO aU ■ ■■ ■■ ■ ■ 

Ci.r'> ' F i't Uti ■ :■ ■ .' 

UIQIUO eWTHSLpy calculation starts USTf^i 

f.i]UAtlO>^ OF CORRESPONDING STATES i:: 


DU 

TCDOt’slCCU 
0Cr'.\3'=UCC i ) 

cpa=cd-/a»aci) 

CC8 = CP'^A?fKI) 
CPC-CPi'APCCi) 
CP5=c?v.\PDa) 


CA!*L Tua'CTCOAT, PCUAT,T,PT, CPA, CPB, CPC, CPD,OHDAT,CPLD ax,. 

I UKDUSFf) 
nfiTUCjSP,I)=f}ELHL 
CPLC JSP,T)^CPLOAT 

TFCK,cy.i:o,'#5 GO to 310 

if}Rltfe:C22^2'^33 JSP,t,BNTLCJSP,l),DeLHVT,CPLCJSP,l) . 

>9 3 F0P4 A t C 5 ■< , X4 , 3X , 'COMP. NO= ' , 1 3 , 3X, 'HL= ' ,El 5 . 8 , 3 X , ' OP VT- ' , f 1 0 . fi , 

1 3X, ^CPL=',K15,8) ■ ■ 

510 CnNTIrJUe . ■ ■ 

500 CO^'TInUE . 

retur^?-j 

E.^0 

subroutine EOUIl.(T,PT,X,Y,KI} 

1 

: variables: 

f STEMP in K 

: , PT; pressure tm 

X,y:MOLE fraction in JilOUIO & VAPOR PHASE , 

: Ki:EaUILIBftIUM RATIOS ' 


. ■ 

; THIS Si,!BPOUTlNE OETEPMINFS 'aiHETHER TO USE UNTFAC , 3RK , PP FUR 

Z FINDING 

: FOUlLl'iPTUN RATIOS AND THEN CALLS THE APPROPRi AXE SUHROIITINE 

SlU^iUliSTlNK EO!JILCT,PT,XXX,YVY,KII) 

OI’*‘-U’Sla.; XiS) , Y(&-) ,KJ CS) ,XXXCS3 ,YYYC6) ,KII(5) 

CO-t U‘Vv f.pAT/ri, ‘,CtjCP,K.SXT{t,ivP,NpPI 
REAL yj , f : i i 
DO lUv/W 1=1 ,N 

xcr)=Kxx(i) ■ ■ . 

1000 Yci)=yYici) 

c 

GO TnC10(.,9,tf:D)KSrTU 
C SRK Af.o P!> PbTHnns 

100 TF CMpn.T.EO.t) GO TO U 

IF (NH.flO, 1 ) GO TO 11 
MRT rr.cPO/^s) 

XKCf.SlTu.r: ), 3) GO TO lU 
aRirFC''!2,3) 

3 F0R'''A r(NX, ,3a, 'BKS-EONS USED FUR CALCULATING EOUlUBRTtJM PAT', 

1 'Tu5',7X,l'^‘S'l 

CO, io in 

ill FOPOAtCS.^! In , 3 X, 'PF.NC; RGUINSON used for calculating eoilxbr 

iiu''! foun Li' ,9<,i'n) 

11 CAH. F')'iTL3C:<,Y,T,PT,KIl,L,KStT!J) 

GO TO L-’ 

9 COUTImiJK 

C IFhC T‘J he used 

1 'f&'? » 4 1 ■ ' ' ' 

4 FnP.'tArC«;Ji,tfl*,3X, 'UfUFAC used for calculating RaUILXBRTUM PATIOS', 

1 ' R ^ , i !'i * ) 

7 Fri5'oA’rC5X,li5^,L3X,lH*/,6X,lri»,63X,lH*/,5X,65CtH*)//) 

10 COTTIO'IK 

Ei-m 




: SUBBOUflNE E0UIl.2CN,f*»COMP,T,PT,J(,y) 

' ¥,’^'■'1 A ■•}•,” 3; 

; iP, ! C’-JTS. FACTION PARAMETERS 

: >;‘;UV:F'J'5ACITy OF the vapor phase 

; ?.iil.:Fij6AC,m Of THE LIOUIO PHASE 


KJ .'■;vlI.lJ^18HIU.« RATIOS 
2;CiiHPRE5SlBIi:.ITy 




THIS PROORA'^ CALCMOATES THE EOWXLISRISjM RATIO 
OR PR EOS’S OEPEMOING UPON KSTTU=10R3 


OS tlSiNG HKS 


2000 


2002 


3000 


3001 

4 

C 

3010 


3020 

3021 


SIIBROOTIMS E0Urij2CX,T,T,PT,KfN,KSITU) 

Di:4Ef^Sl'a<^ X(5),¥C5),K1C5{5) 

P 1 4EMSin.‘4 A b K K c 5 5 I PHt V ( s 5 , PHIL C 5 ) , XX I C s ) , YY I ( 5 ) 

COMMO^/&YSOAT/PCC5),TC(53,ZC(5),OMeG4(5).VCCS) 
COMMON/EQIDAT/F C 5 ) , B ( S) , FM , SM , A. A C 5 ) , AO , \htMX , BS T A R 


KICX r JlJlHTERACTION PARAMETERS 

DO Pll I»lr« 

DO 9U^J5J,« 

ih^Mur-^ 


SoilfSHEM THIS SimROUTINE IS EXECUTED FOp VAPOP 
Sfitmsn Tars SOBROOTTf-E is executed for liouio 
HOaO 
ReS2.05 

SU0ROUTIHE ROUT CAOOED CHCE FUP 010. Ar^O OnCF- EOS VAOOP 
C A I. L R 0 0 T C Z .«! A K , Z Ml fl , E 1 1 Y , T , P T , z , H , n a , K s I T iJ 1 
WRirFC22,H31) ZMAX,ZOiM , , ,, 

Fn»V.Ar(5X/ *ZMAX=' .F7.4, 'ZMIHs:%F7,4) 

IF C CZmAX,LT.O.O) ) GUTO 9990 
IF CD'1. .E.l 3 GOTO 12 

FOR LIUUTO PHASE ZHIM CORRESPONDS TO Z-fAX 

ZiiAXsZOTT 

VsYOfilMG 

V=CC?,-.«X)^-R*T/PT3 

CD TO (2U'->«£iOOn,20'}23K.SlTU 

GA=R, 1274703 

GB=C,nRGbl03 

DO 2 r = l ,-4 

A(T) = CGA»H*T»riCI)*F(X)/GB) 

RM=CGAl=P’»'T*BM^Fn/GB3 

nn 3 1=1, ' ■ . 

RAHii = 0,,: 

DO 4 

GO 10C3«.'R'>3, I, J<'P3 )KSIT!1 

HABA=linT,AtC(f ..3“1‘HI,J3 3*U A(X)*A(J) 3**0.5>»YCJ)/AM) 

GO ID 4 

HAS A=H A B A +2 , * 1 1 . -Kt ( J , ! 3 ) *SORT ( A A C J ) *AA C I 3 ) * Y C J 3 / AM 

COi^t'lMUE , . ' 


PrilLfi:r..ionti» fugactty 

GO iC! C'J'hu.lOOC , 3U203KSTT0_ 




p H I L N = C H C L ) ^ '( Z M A i - 1 ,' 0 3 / H M 3 + A D I'! G C V / C 2 M A X ♦ C V - H M ) ) ) + c A H / C 8 : 1 * B * 1' ) ) 
l«Al,DGC(V+B,-U/V)^(CHCI)/BM3-2.0*BABA) 

Ff’ TD 

PHlU4=nC i' ) *'(ZF/\X-1.3/BM-ADOCCZMAX-0STAR)-ASTAR/C2.a28#BSfAP3 
l^CdABA~HCl'3/H.M3*ALkcCZMAX + 2:414*BSTAP3/(ZMAX-0.4H*BSTAR3 3 

cod IT 'I OF. 

IF tUO.EO. 1 ) GOTO 5 

0£PEh*)TM.: 3 Df>i ?nD ; ASSIGN THE FUGACITY 
PHIVC I 3=EXPCPHILU) 

GOTO i 
COK'TIUUF 

pHlUCDxEXPCPHrLM) 

COUTIi'OJE 

ISipiAT^tHF. DUOP FOR THE LIQUID 

yHW ... 

Guro to , : ' 






CAr4C'ffi;n'in'’‘4 OF BOTH LIQUID AND VAPOR fUGACITY COMPLfcrrF 
x, i 

DO 1 R 

Y ( I j = X T f/t ) 

Yrn=y/i(!) 

COi'fl'T'ifji-; 

[).-) I S f. = 1 ,?■] - , 

if' U:F(2:i, i'/it) I, PHIL( I) ,I,PHIVCI) 


^ ** f » A. 4 I - \ 4 . r A ' a i ^ A $ r n 4 w \ 4 / # * . r f 1 * ¥ / . . 

1731 FUB:iAT(5X,'PaTLI0(''.I2,^)s',Ei6.9,'PHlVAP(»,T2r') = ‘»F7,5/} 

fC(l) = (PBTLCl)/PHIV(i)3 i V , ^ 

•.?RiT^;c22,ta2),T.Ka) 

102 FO^^'IATC' KCM2, ^J5%F15,93 

11 cori loiio 

GO'"';) 1000 

If i:‘H0 mM> ROOT XSi WE CASE OF OME PKAL IS MKOATIVO -PPI.:;T 
: ERROR S-SIOP 

?999 .^RTTFC22,35) 

>5 FOR-'iATC' ERftORiP#$%"& CHECK?? CHECK ????? '5 

STOP 

iOOO COr^ITNUF 

RFTUP.1 
F;jO 

! SUHROUTI«E ROOT(Zf4AX,ZMIH,Kl ,I,T,P,X,r<,NO,'<.sn‘in 

: VARIABLES: 

: tr:hel>ijceo temperature 

: PR: REDUCED PRESSURF 

: THIS SUBROUTTRE calculates the compressibility FAC7UH,7, by 

FVALUATIMG the routs of the cubic rso.JATIuO AND i'HEM ASSIGHTfiG 
: the largest real root to the vapor Phase and the smallest 

: ROOT TO TrtB LIQUID PHASE 


subroutine RnnT(7.MAX,ZMrR,Kl ,y ,T,P/Z,L,flOYKSIT!U 
oirtEMSins y(5),TPC5) ,PH(S) 

DIMFNSIOH TFMPCS) ,icuy,5) 

C0HM0N/5YSDAT/PCCS) .TC(S) ,ZC(S) fOMFOACS) ,vrtS} 

CO.MMD •: /E Q n AT/F C 5 ) , cu 5 ) , FH , B?-’ , A A ( 5 ) , \M , ASX AR , tiS TAR 

PEAL K 1 ' 

R:GAS rnr3r-\?;i:s 

p ^ 0 I , 

IF (ysiTu, v'V.II GO T{) 2000 


SRK EQ'LS JStB, 


iITU = t 




G A f v’S • C ',t .1 S 5’ t\j‘'* fS 1. M RKS 

GA='y 42 7 

FMt=V,4< 

CALC Hf. ATE RFOUCEn rFOPEPATURF & PRESSURE 

PFL = t .(• ■ 

Df’ 1 1 = 1 , ■) 

TpCi)=T/rC(1 ) 

PR (n=p/p: Cl. ) 

CAtX'Ur.ATT LlFr5:RFfJ'T VARIABLES OF RKS 

B(I)=GB’fP*’TC(l)/(PC(X)3 

B..^ = B'^ + Y(T)I'HCi) 

TKHPCl3 = C5 ,''' + C0,lBU + ! ,5744nMFGAtI)-0.l7b*0MFGACl)*0MEGA( Ul-^Ci ,0 

1-C (TRC 1 ) ) ) 3 

FCT) = (l,s#/TBCl))l((TKMPCl))’^*2.0) 

CruTnjOif.; 

CALCtn.AlTOL' OF VAPJJ^BLES (COM.) 

DO 2 1 = 1 ,1 . ' 

on J ii s 1 - 

rfnl=Y(.1 )*V CG)^( tCtCCD^FC I)^=TC(J)*F(J))/CPCtI3*PCtJ) ) 
tO-K! ( J,J 3 )+FMj 

cooriLt.iL 

DEFSOEX+CCYCI) yfCCl3l/PCC.l}) 

COLTIluIK 

Calculh nnu of vapiahles ccout.) 

^HsCpoi/ncN) 
eSTAR=CBH<^P3/(P.*T3 . 

ASTAPsCC-LSTMRtC.M/'tlO+FO 

CALCULArTLG RiiutS OF THE CUBIC EQUATION 

H 1 ss-** I f V ,■'■■■■'"■■■■ ■"■ ' ' ' ' 

ClsCRSi’AH-Ha'IAK-CBSTAR^flSTAR)) 

Dl»*c ASTAR^I'OSTAR) 

GO TO 2m, I ■ ■ ■ 

Coutinhe 

■ ' ■ ■ ■ . ■;■■■ 

DC auO? ( = !,« ', ■ ■ 


2002 

2003 


200i 


155 

iO 

c 

r 

^5 


C 

c 


111 


c 

12 


c 

20 

c 

21 


'&66 

667 


. C 

t c 


>9 



B C 1 3 = f ; , ■') 7 7 3 0 * p r C c 1 3 / PC C I } 
e '5=L-B’HfUi3nrci3 
HAf'i )=sfc, **572ni 
AI.(I)«C (1 ,nu.3, .. 
n=S‘tl,“((TftCI)3^>*= 

Crc^Ti.PlE 

DC 2*^0 3 1*1 , N 

BKTAHsari*P/(ft*T) ' 

Rt=BSrAR-|, ,0 

C 1 = A S T A H - 3 . 0 C 8 S T AR ♦ aSf AR 3 - 2 . 0 =»= BST A R 

DI=-RSTAR*-CASTAR-BSTAR-(BSTAR*fBsfAR3 3 

CU^XT:v'i,}(;; 

PP*f3,0*Cl7CBl *81 33/3,0 
!30-(2/,O*ai-9,0*Bl*Cl+2*O*81*Bl*B! )/27. 
feRs(Pplpp#pp3/27.0+Q0*05/4.0 
RpIDeTKR^SpF:S WHICH PATH t 6 TAKK 
IP {CA8SCnR3).GT.O.iB*203 GOTO 155 
(j Cl X 0 2 0 

IF (RR) 10,20,30 
CrWTlM&C 
WRITS (22, 95 3 
R<0 

PORMATC* THERE ARE THREE UNF;0UA6 REAL ROOTS OF TH.: 

i u s*) 

THIS GIVES IM RADIANS 
PHI»AC0SUS3 

AMP®2,y*a'’PP/3.0)**0.5) 

A«G=(?HI/3,0) 

UNPOUAE ROOTS 


0,26992*OHECA ( I ) *iV^.F/,;n ( r ) 


0 


CI 58 TC cMD , 



^0943953 
:3«AMP*Cr3v5(ANG+4. 18879 3 ' 

IF (OQ.OT.O.O) GOTO Ut 

XlsXi”CB3 /3.0) 

X2=X2-{Bi/3.03 
X3 = X3-CBl/3,<)) 

»'PITFC22,12) X1,X2,X3 , . ' 

for AT C ' PDOTR ARE : * , 1 OX , ' X 1= ' , FI 0 . 6 , :3X , ' X2= ' , Fl 0 . 6 , 3 X , ' X3= ' , F 1 0 

Ff«F5 XtiF :AX1MUM & B.JtijMllK OF READ HOuTS 

CALI) Z'/Ar.'.’^CX!. ,X2,X3,’if'AX,7„'V!j!0 

E -= T‘".AX 

WnirF(22,4R),Z 

GOTO n 

for p=m 

CC-jXI.hUI', . 

WRITE (22 ,21. 3 

Fv'iPB.AT(' THFRt*, ftPE tiiPEf, HEAL F0f)fKC2 FOlJAL) OF THE CUBIC FOT.i 
1 , 1 D * 3 

IF (0/.GT.0,0) GOTO un6 

AX=-( CO..i/2.i.'3'ff (1 ,n/3,U) ) 

GOTO Ph7 

AX-C(30/3.u)**Cl ,^'/3,rs)) 

Xi*(MX + AX)-CPS /3.'*3 

X3 = ’-AX-(iU/3.r) 

X2=X3 

Fi^OiriC r.fE MAX, * .V.I1, OF DIFFERENT REAl. ROOTS 
CALb /tfhL.:'l::(:a.,>'.2,X3,^-lAX,E llw) 

Zs:Z/AX., 


•. ii .* TTF :(2/., r /), Xl ,;(2 ,Xi 
WRITE' (22, 49) ,Z 


CCTO I 1, 

P>0 
CO^TlbHF 

WRirP(2? ,313 

FORhACC' THi; cubic EOUATIoh 1h Z HAS ORE REAL ROOT TWtl COMPLEX 
I ROOTS ') 

CAbCUlU'cnOM OF THE Ch-'H:, heal root ■ 

AX=(-CU0/2.O3'HSUFTCftR)l 
IF (M,GT,0,v) OTTO 668 
AXa-AX 

AXa-CTX#f{jl.O/3.a3) 

GOTO 669 ^ ■■ . 

AX« ( Aa** C 1 ,0/3 .0) ) , • 

I|*"c|x*».nu0/ 3.C3 3 . 

finTii H7 5 t ^ 


3C01 


3002 
300 3 

1621 FOR“-lATC5X,"r{If'!CTHl.'i=',f.if',0) 
fi£fOP»l 

EIj’O 

c 

C ■. ,,. RlJSHOtirXWE Z¥Aftiff:(XI fX?rX 3 .Z«AX://,’^I^I) 

c 

c 

C fHIS SUBROUTINE CALC'J-LATES THE MAXIMU.-I ANl.) , MINIMUM OF THREE 

C . ' IW no. XI ,X2,X3 

,, '' '■ ’iuBfiOOfINE ZVALUF.CXi ,X2,X3,ZMAX,ZMIisi)' ^ ' . ' 

r khL THREF POSITIVE 

IF C(X1 .GT.d,0).A:^!).CX2,GT,0,0).AND.CX3.0T,0,0)) GO TO 10/ 

C ML TdPES ^JKGATIVE , , ^ ^ ■ 

If (fxt ,1,1, 0.03. A GJ,(X2.LT. 0,0). Af^D.(X3,l.T, 0.0)3 GO TO 20 

C IF ]>'□ A3G POSITIVE 

If C (XI .liT.O.03 GO TO 30 

C (:ASE4:lPnari'i:VCe2 i‘.EGATIVE bouts 

ZViAX = AHftXl (XUX2/X33 

7f.•I:v=7,"^X 
GO 13 S 

C CAS'-t:HUL POSITlVii 


BX=:CHX*>i-’( I ,0/3.0)) 

ZsAA + iiX-’C 41 /3,P) 

IF C'-, Cl. GO ru iso 
ZHnA=z 

THF Pc'.'.Ji iDilv IS ASSCGO^CD TO ZHAX nn Oiv. 

.. ARC fJuATI.,1/ IS s)(;i.!'ih! Fnp VAP!p-( OH TOE ^jTjHio v,,k 

a OTPnTHKT1.C/iO VAUU£ Ao fHtS IS OO-T 

GIVE VAUOK HYPOTHKTICAU-IO 

2MT .(' 

CO 'I’O ?0:C 

Cr/OnriOr 

ZHfOsZ 

GIVE HAXI'/n>5 VAflijE HYPCTHETICAU 10 

cn'^ri’-KfE 

/*F,'IIFC22,49),Z 

FCP^IATC' oof: real BOOT=',MO.n 
C0WTI«U>; 

FOPfiAlC*^ ‘THE VALUE OF Z TAKCK :',F!5.p) 

GO TO (300 1 , 300 1 , 300 2 ) BS ITU 

F I M = C z * * 3 . y 3 - / ♦ 'i + Z * c A s :r - s S 1 A P - B S 4 ' ■ * B .5 T A n - T S r A k t y. S T -'i H 
GO i'3 3 OP 3 

Ff«sCz#% 3 . 0 )+Fl*' 2 *Z+CHZ + I>l 


10 

Z?«AXsfv'''AXl(>.l.X2,X31 
Z''iTOsH-r..'l (X3 ,X/,X3) 
GO TO 5 


c 

TASK?; AS, L ■■'FGATlVr; 


20 

xrt n:i7^,^^) 

Fr!P''?AT ( s / , ' F.RRnR-ALL 

GO ro 5 


7 

z nfgattvf') 

C 

Ci\SK3: A’.rJi'i U.'iSfTlVK S 

OX'E HAGATIUF 

30 

IF’ ((>;•! -K'/)) .p) 

ElsXI 

A2=X? , 

GO TD si; 

;!j TP -IP 

40 

IF (CKl'+Xi). { 
Al=A‘i 

A2FX3 

G!1 T') 50 

if* TC' 6 0 

60 

A 1 s X 2 

A2 = X? 


50 

ZMAXS'I-UAXI (At ,h2) 


5 

CClPTMUC 



Rf;T!JR?( 



r ^ $ 3 ^ ^ ^ f * f 4 f f f* # # f f f -t * f t -f * f * f ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

SUBH-.IUTI’-'. /. 0 ;/T( 7 / 1 AX//T'\TTi,Y,F,P jM>'!,rC,PC,TfPiU*^El 34 f'X,. 1 ,KSlTt.n 
DlMaMSCOi/' Y(25 3 , ?C C ’>‘> 3 , PC ( 2^3 ,Q<1 EGa (25) ,tR(2h) ,P(H2:)3 ,B(253 
DIME f>( S T ft A T r: ''4 P C 2 5 ) , F ( 2 5 J 

C #as sjVF/tui’rcE giJks z GJVEP T.p,v usi^c sax na fr fur ctae 

e , pASErVAPOiin ALn;^;. Z IS THE CORRlkt COPRSSSTBILTj: FACrnR 

f . ' §2^5ptZ,i i.'/. It ' r. .* t "I A 1 


1001 


1002 


1003 


1004 


■r; A“‘ ' • - / 7 4 s/ 'i ■■ ■ ■ ■ ' ■ ■ ■ 

r,’-!=*' .<•! ■ ■■ ■ ■ ■ 

F '.1 

OEMS'"’. 0 

DO 1^:1 A! . . ■■ .A 

TP(I)=T/tCCI) 

)s'voncn ■■■ '■ 

B(ij=',jB’»'R*Tccr)/(pccn) 

Brt=M.vt+ycT)i-Ha) 

i)*n,.iKu-^( C j vm J 

1 - ( U H Cl ) ) ^’*0 ,5 ) ) > 

F(TJ*{l.0/1’HCT))*CCfEMPCI))*»2.O) . 

C'FiiifjME : 

on 2 1 = 1 , ?| 

DO 3 J=1,M 

f Mi=v(i)=»‘y(jHcc CTCct)»r(i)’!'TC(a)»FCJj)/cpcci)»pc( j) 1 )+?■:- 

COOTFf'HJt' 

DF:?lsDg?A' + (CY(I3>f!TCCI>)/PCf D) 




BSTAH=CHM»F3/tR*T3 

ftSTAR=at^C>TiR*C3A5/GB)*FM 

ClsCArjTAP-'BSTAR-(0STAB»8STAR)3 

Dl=-CAr>TA!3*BSTAR) 

Go ftl I0u4 

BM»0, 

FM=0, 

00 1002 Is! 

TR C I ) sf /IC ( 1 ) 

B(I)*0.O778O^R*TC(I)/PCCI3 
RHSBM + BCDfYCI.) 

ALPHAl = l,f CO.374644-1 -542?b*C)MEGA{ n-0.2t3‘»P2*0-1fc-GA( ! )f')MKGA( 1) 

i)*ci,-C(TRcn)*^o.55 ) 


hhPHklsiiWtlkl^khPnM 
r(I2»C0.4S724*C(R^TCCX5)^*2)/PC{I}3^AGPHAI 
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00 :100| 

00 i®^03 JSI.N 

PMsF‘«+Y(I3*f {J3*((FCr)*F(J))*»0.5) 

ASTAR=FM'^'P/(Btt*P>T3 

RSTAHsBM4:p/CRtt) 

BlsfiStAB-UC 

Cl=:ASTAR-3.n*CBSTAR*BSTAP3-2,P*B5TAR 
Dls«B.STAK*(ASTfiR-BSTAR-(PoTAR’^BSTAH) 3 

PP = p.wfCl-('>!*bl ))/3.0 ^ ' ■ 

00s(2? ,n*i) 1-9.0*01 !t«Cl +2. 0»B I *E 1*813/27.0 
rjps(Pp*,--t.>*t»o)/27 ,0 + 00*00/*,') 

If ( (.\ ^i.(t.M3).Gl'.''.lF.-103 GOTO 155 . ■ 

ppsf.c 

IF CEP.) 1..>,20,3C' 

COEIMmU!,:. 

FORPAfC' T-MtRK AhF THnecl t/MEOEAL ROOTS OF ThE CUBIC F^T'^. T^( /' 
XS=CCO ) + 0, 1/1. :r)/CC-( ?P*?P*PR 3 3/27,0)) 

KSsSORTCXS) 

F + 3=).COSCXS) 

.AwP = 2.0*C C»OH/3.O)**0.5) 

ANGs(0-ti/3.0) 

XisAMR* (C-jSC/ifiG) 3 , 

K2 = A'’P*rnSC A.'it.; + 2.u0 4?95) 

XJsAMp^cOS C i’iriG + 4 .1 8fe?9) 

£F COO.ET.O.O) GOTO 1.31 
Xl=-/l 
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1*5) 
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Xll|!!l?fl3'o^527 
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TERMIsRl'ALOGC CV-3)/V3 

rERM|sCA/CB*Si)RTCT) ) )^ALOG( (V+P)/7) 

tERfi|3«R*ALnGCV/70) 

OELA*»CfeRMt»T)~fERM2-CT*TKRM3) 
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CPL = CP')f CVXPPS'rt j 
fc:P ■*? = (! 3^4 

DELnv'7=C( ?.(.‘6»TF?:--3)3ClC.P5't-'nMFCA*TFR.M4))*R*TC 

f*! 0 
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